max/poet
1
0
Fork 0
mirror of https://git.gfz-potsdam.de/naaice/poet.git synced 2025-12-15 12:28:22 +01:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'ml/new-iphreeqc-api' into 'main'

Browse Source 
Update to IPhreeqc v3.8.2

See merge request naaice/poet!37
This commit is contained in:
Max Lübke 2024-11-07 14:23:58 +01:00
commit 66098aab40
8 changed files with 292 additions and 224 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
R_lib/init_r_lib.R
View File

@ -13,6 +13,34 @@
### this program; if not, write to the Free Software Foundation, Inc., 51 ### this program; if not, write to the Free Software Foundation, Inc., 51
### Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. ### 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 pqc_mat matrix, containing IDs and PHREEQC outputs
##' @param grid matrix, zonation referring to pqc_mat$ID ##' @param grid matrix, zonation referring to pqc_mat$ID
##' @return a data.frame ##' @return a data.frame
@ -33,7 +61,7 @@ pqc_to_grid <- function(pqc_mat, grid) {
res_df <- as.data.frame(result_mat) res_df <- as.data.frame(result_mat)
# Remove all columns which only contain NaN # 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 # Remove row names
rownames(res_df) <- NULL rownames(res_df) <- NULL

2
ext/iphreeqc

@ -1 +1 @@
Subproject commit e6e5e0d5156c093241a53e6ce074ef346d64ae26 Subproject commit 61214a01ad4cf99527f657e6a41c68282e4886c4

9
src/Chemistry/ChemistryModule.cpp
View File

@ -1,6 +1,7 @@
#include "ChemistryModule.hpp" #include "ChemistryModule.hpp"
#include "PhreeqcEngine.hpp" #include "PhreeqcEngine.hpp"
#include "PhreeqcMatrix.hpp"
#include "SurrogateModels/DHT_Wrapper.hpp" #include "SurrogateModels/DHT_Wrapper.hpp"
#include "SurrogateModels/Interpolation.hpp" #include "SurrogateModels/Interpolation.hpp"
@ -167,9 +168,11 @@ poet::ChemistryModule::ChemistryModule(
this->n_cells = chem_params.total_grid_cells; this->n_cells = chem_params.total_grid_cells;
if (!is_master) { if (!is_master) {
for (const auto &[pqc_id, pqc_config] : chem_params.pqc_config) { PhreeqcMatrix pqc_mat =
this->phreeqc_instances[pqc_id] = PhreeqcMatrix(chem_params.database, chem_params.pqc_script);
std::make_unique<PhreeqcEngine>(pqc_config); for (const auto &cell_id : chem_params.pqc_ids) {
this->phreeqc_instances[cell_id] =
std::make_unique<PhreeqcEngine>(pqc_mat, cell_id);
} }
// for (std::size_t i = 0; i < chem_params.pqc_ids.size(); i++) { // for (std::size_t i = 0; i < chem_params.pqc_ids.size(); i++) {
// this->phreeqc_instances[chem_params.pqc_ids[i]] = // this->phreeqc_instances[chem_params.pqc_ids[i]] =

25
src/Init/ChemistryInit.cpp
View File

@ -2,16 +2,14 @@
#include <Rcpp.h> #include <Rcpp.h>
#include <cstddef> #include <cstddef>
#include <fstream>
#include <string>
#include <utility>
#include <vector> #include <vector>
namespace poet { namespace poet {
void InitialList::initChemistry(const Rcpp::List &chem) { void InitialList::initChemistry(const Rcpp::List &chem) {
this->pqc_solutions = std::vector<std::string>(
this->transport_names.begin() + 3, this->transport_names.end());
this->pqc_solution_primaries = this->phreeqc->getSolutionPrimaries();
if (chem.containsElementNamed("dht_species")) { if (chem.containsElementNamed("dht_species")) {
this->dht_species = Rcpp::as<NamedVector<uint32_t>>(chem["dht_species"]); this->dht_species = Rcpp::as<NamedVector<uint32_t>>(chem["dht_species"]);
} }
@ -69,21 +67,14 @@ InitialList::ChemistryInit InitialList::getChemistryInit() const {
// chem_init.field_header = this->field_header; // chem_init.field_header = this->field_header;
chem_init.database = database; 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_scripts = pqc_scripts;
// chem_init.pqc_ids = pqc_ids; // chem_init.pqc_ids = pqc_ids;
for (std::size_t i = 0; i < pqc_scripts.size(); i++) { // for (std::size_t i = 0; i < pqc_ids.size(); ++i) {
POETInitCell cell = { // chem_init.pqc_input[pqc_ids[i]] = pqc_scripts[i];
pqc_solutions, // }
pqc_solution_primaries,
Rcpp::as<std::vector<std::string>>(pqc_exchanger[i]),
Rcpp::as<std::vector<std::string>>(pqc_kinetics[i]),
Rcpp::as<std::vector<std::string>>(pqc_equilibrium[i]),
Rcpp::as<std::vector<std::string>>(pqc_surface_comps[i]),
Rcpp::as<std::vector<std::string>>(pqc_surface_charges[i])};
chem_init.pqc_config[pqc_ids[i]] = {database, pqc_scripts[i], cell};
}
// chem_init.pqc_sol_order = pqc_solutions; // chem_init.pqc_sol_order = pqc_solutions;

160
src/Init/DiffusionInit.cpp
View File

@ -1,4 +1,3 @@
#include <algorithm>
#include <tug/Boundary.hpp> #include <tug/Boundary.hpp>
// leave above Rcpp includes, as eigen seem to have problems with a preceding // leave above Rcpp includes, as eigen seem to have problems with a preceding
// Rcpp include // Rcpp include
@ -8,7 +7,7 @@
#include "DataStructures/Field.hpp" #include "DataStructures/Field.hpp"
#include "InitialList.hpp" #include "InitialList.hpp"
#include "PhreeqcInit.hpp" #include "PhreeqcMatrix.hpp"
#include <Rcpp/Function.h> #include <Rcpp/Function.h>
#include <Rcpp/proxy/ProtectedProxy.h> #include <Rcpp/proxy/ProtectedProxy.h>
@ -53,98 +52,102 @@ static std::vector<TugType> colMajToRowMaj(const Rcpp::NumericVector &vec,
} }
} }
static std::vector<std::string> // static std::vector<std::string>
extend_transport_names(std::unique_ptr<PhreeqcInit> &phreeqc, // extend_transport_names(std::unique_ptr<PhreeqcInit> &phreeqc,
const Rcpp::List &boundaries_list, // const Rcpp::List &boundaries_list,
const Rcpp::List &inner_boundaries, // const Rcpp::List &inner_boundaries,
const std::vector<std::string> &old_trans_names) { // const std::vector<std::string> &old_trans_names) {
std::vector<std::string> transport_names = old_trans_names; // std::vector<std::string> transport_names = old_trans_names;
std::set<int> constant_pqc_ids; // std::set<int> constant_pqc_ids;
for (const auto &side : tug_side_mapping) { // for (const auto &side : tug_side_mapping) {
if (!boundaries_list.containsElementNamed(side.second.c_str())) { // if (!boundaries_list.containsElementNamed(side.second.c_str())) {
continue; // continue;
} // }
Rcpp::List mapping = boundaries_list[side.second]; // Rcpp::List mapping = boundaries_list[side.second];
const Rcpp::NumericVector cells = mapping["cell"]; // const Rcpp::NumericVector cells = mapping["cell"];
const Rcpp::NumericVector values = mapping["sol_id"]; // const Rcpp::NumericVector values = mapping["sol_id"];
const Rcpp::CharacterVector type_str = mapping["type"]; // const Rcpp::CharacterVector type_str = mapping["type"];
if (cells.size() != values.size()) { // if (cells.size() != values.size()) {
throw std::runtime_error("Boundary [" + side.second + // throw std::runtime_error("Boundary [" + side.second +
"] cells and values are not the same " // "] cells and values are not the same "
"length"); // "length");
} // }
for (auto i = 0; i < cells.size(); i++) { // for (auto i = 0; i < cells.size(); i++) {
if (type_str[i] == "constant") { // if (type_str[i] == "constant") {
constant_pqc_ids.insert(static_cast<std::size_t>(values[i])); // constant_pqc_ids.insert(static_cast<std::size_t>(values[i]));
} // }
} // }
} // }
if (inner_boundaries.size() > 0) { // if (inner_boundaries.size() > 0) {
const Rcpp::NumericVector values = inner_boundaries["sol_id"]; // const Rcpp::NumericVector values = inner_boundaries["sol_id"];
for (auto i = 0; i < values.size(); i++) { // for (auto i = 0; i < values.size(); i++) {
constant_pqc_ids.insert(static_cast<std::size_t>(values[i])); // constant_pqc_ids.insert(static_cast<std::size_t>(values[i]));
} // }
} // }
if (!constant_pqc_ids.empty()) { // if (!constant_pqc_ids.empty()) {
constexpr std::size_t keep_h_o_charge = 3; // constexpr std::size_t keep_h_o_charge = 3;
for (const auto &pqc_id : constant_pqc_ids) { // for (const auto &pqc_id : constant_pqc_ids) {
const auto solution_names = phreeqc->getSolutionNames(pqc_id); // const auto solution_names = phreeqc->getSolutionNames(pqc_id);
// add those strings which are not already in the transport_names // // add those strings which are not already in the transport_names
for (const auto &name : solution_names) { // for (const auto &name : solution_names) {
if (std::find(transport_names.begin(), transport_names.end(), name) == // if (std::find(transport_names.begin(), transport_names.end(), name)
transport_names.end()) { // ==
auto position = // transport_names.end()) {
std::lower_bound(transport_names.begin() + keep_h_o_charge, // auto position =
transport_names.end(), name); // 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 // static Rcpp::List
extend_initial_grid(const Rcpp::List &initial_grid, // extend_initial_grid(const Rcpp::List &initial_grid,
const std::vector<std::string> &transport_names) { // const std::vector<std::string> &transport_names) {
// std::vector<std::string> names_to_add(transport_names.begin() + old_size, // // std::vector<std::string> names_to_add(transport_names.begin() +
// transport_names.end()); // 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<Rcpp::List, Rcpp::List> std::pair<Rcpp::List, Rcpp::List>
InitialList::resolveBoundaries(const Rcpp::List &boundaries_list, 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 bound_list;
Rcpp::List inner_bound; Rcpp::List inner_bound;
Rcpp::Function resolve_R("resolve_pqc_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->transport_names = extend_transport_names(
this->phreeqc, boundaries_list, inner_boundaries, this->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) { // if (old_transport_size != new_transport_size) {
this->initial_grid = // this->initial_grid =
extend_initial_grid(this->initial_grid, this->transport_names); // extend_initial_grid(this->initial_grid, this->transport_names);
} // }
for (const auto &species_name : this->transport_names) { for (const auto &species_name : this->transport_names) {
Rcpp::List spec_list; Rcpp::List spec_list;
@ -179,8 +182,11 @@ InitialList::resolveBoundaries(const Rcpp::List &boundaries_list,
c_type[c_id] = tug::BC_TYPE_CLOSED; c_type[c_id] = tug::BC_TYPE_CLOSED;
} else if (type_str[i] == "constant") { } else if (type_str[i] == "constant") {
c_type[c_id] = tug::BC_TYPE_CONSTANT; c_type[c_id] = tug::BC_TYPE_CONSTANT;
c_value[c_id] = Rcpp::as<TugType>( c_value[c_id] =
resolve_R(this->phreeqc_mat, Rcpp::wrap(species_name), values[i])); static_cast<TugType>(phreeqc_mat(values[i], species_name));
// c_value[c_id] = Rcpp::as<TugType>(
// resolve_R(this->phreeqc_mat, Rcpp::wrap(species_name),
// values[i]));
} else { } else {
throw std::runtime_error("Unknown boundary type"); 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++) { for (std::size_t i = 0; i < inner_row_vec.size(); i++) {
rows.push_back(inner_row_vec[i] - 1); rows.push_back(inner_row_vec[i] - 1);
cols.push_back(inner_col_vec[i] - 1); cols.push_back(inner_col_vec[i] - 1);
c_value.push_back(Rcpp::as<TugType>(resolve_R( c_value.push_back(static_cast<TugType>(
this->phreeqc_mat, Rcpp::wrap(species_name), inner_pqc_id_vec[i]))); phreeqc_mat(inner_pqc_id_vec[i], species_name)));
// c_value.push_back(Rcpp::as<TugType>(resolve_R(
// this->phreeqc_mat, Rcpp::wrap(species_name),
// inner_pqc_id_vec[i])));
} }
inner_bound[species_name] = inner_bound[species_name] =
@ -276,7 +285,8 @@ static Rcpp::List parseAlphas(const SEXP &input,
return out_list; 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 boundaries;
Rcpp::List inner_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)]; diffusion_input[DIFFU_MEMBER_STR(DiffusionMembers::ALPHA_Y)];
const auto resolved_boundaries = const auto resolved_boundaries =
resolveBoundaries(boundaries, inner_boundaries); resolveBoundaries(boundaries, inner_boundaries, phreeqc);
this->boundaries = resolved_boundaries.first; this->boundaries = resolved_boundaries.first;
this->inner_boundaries = resolved_boundaries.second; this->inner_boundaries = resolved_boundaries.second;

167
src/Init/GridInit.cpp
View File

@ -1,11 +1,13 @@
#include "InitialList.hpp" #include "InitialList.hpp"
#include "PhreeqcInit.hpp" #include "PhreeqcMatrix.hpp"
#include <RInside.h> #include <RInside.h>
#include <Rcpp/Function.h> #include <Rcpp/Function.h>
#include <Rcpp/vector/Matrix.h>
#include <Rcpp/vector/instantiation.h> #include <Rcpp/vector/instantiation.h>
#include <cstdint> #include <cstdint>
#include <fstream>
#include <map> #include <map>
#include <memory> #include <memory>
#include <regex> #include <regex>
@ -15,37 +17,41 @@
namespace poet { namespace poet {
static Rcpp::NumericMatrix // static Rcpp::NumericMatrix pqcMatToR(const PhreeqcMatrix &phreeqc, RInside
pqcScriptToGrid(std::unique_ptr<PhreeqcInit> &phreeqc, RInside &R) { // &R) {
PhreeqcInit::PhreeqcMat phreeqc_mat = phreeqc->getPhreeqcMat();
// add "id" to the front of the column names // PhreeqcMatrix::STLExport phreeqc_mat = phreeqc.get();
const std::size_t ncols = phreeqc_mat.names.size(); // // PhreeqcInit::PhreeqcMat phreeqc_mat = phreeqc->getPhreeqcMat();
const std::size_t nrows = phreeqc_mat.values.size();
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++) { // phreeqc_mat.names.insert(phreeqc_mat.names.begin(), "ID");
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); // 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, // Rcpp::colnames(mat) = Rcpp::wrap(phreeqc_mat.names);
const Rcpp::NumericMatrix &grid) {
Rcpp::Function pqc_to_grid_R("pqc_to_grid");
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<int, std::string> static inline std::map<int, std::string>
replaceRawKeywordIDs(std::map<int, std::string> raws) { replaceRawKeywordIDs(std::map<int, std::string> raws) {
@ -60,24 +66,25 @@ replaceRawKeywordIDs(std::map<int, std::string> raws) {
return raws; return raws;
} }
static inline uint32_t getSolutionCount(std::unique_ptr<PhreeqcInit> &phreeqc, // static inline uint32_t getSolutionCount(std::unique_ptr<PhreeqcInit>
const Rcpp::List &initial_grid) { // &phreeqc,
PhreeqcInit::ModulesArray mod_array; // const Rcpp::List &initial_grid) {
Rcpp::Function unique_R("unique"); // PhreeqcInit::ModulesArray mod_array;
// Rcpp::Function unique_R("unique");
std::vector<int> row_ids = // std::vector<int> row_ids =
Rcpp::as<std::vector<int>>(unique_R(initial_grid["ID"])); // Rcpp::as<std::vector<int>>(unique_R(initial_grid["ID"]));
// std::vector<std::uint32_t> sizes_vec(sizes.begin(), sizes.end()); // // std::vector<std::uint32_t> sizes_vec(sizes.begin(), sizes.end());
// Rcpp::Function modify_sizes("modify_module_sizes"); // // Rcpp::Function modify_sizes("modify_module_sizes");
// sizes_vec = Rcpp::as<std::vector<std::uint32_t>>( // // sizes_vec = Rcpp::as<std::vector<std::uint32_t>>(
// modify_sizes(sizes_vec, phreeqc_mat, initial_grid)); // // 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) { static std::string readFile(const std::string &path) {
std::string string_rpath(PATH_MAX, '\0'); std::string string_rpath(PATH_MAX, '\0');
@ -98,10 +105,26 @@ static std::string readFile(const std::string &path) {
return buffer.str(); return buffer.str();
} }
void InitialList::prepareGrid(const Rcpp::List &grid_input) { static Rcpp::List expandGrid(const PhreeqcMatrix &pqc_mat,
// parse input values const std::vector<int> unique_ids,
Rcpp::Function unique_R("unique"); 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 script;
std::string database; std::string database;
@ -126,8 +149,9 @@ void InitialList::prepareGrid(const Rcpp::List &grid_input) {
} }
this->database = database; this->database = database;
this->pqc_script = script;
Rcpp::NumericMatrix grid_def = Rcpp::IntegerMatrix grid_def =
grid_input[GRID_MEMBER_STR(GridMembers::GRID_DEF)]; grid_input[GRID_MEMBER_STR(GridMembers::GRID_DEF)];
Rcpp::NumericVector grid_size = Rcpp::NumericVector grid_size =
grid_input[GRID_MEMBER_STR(GridMembers::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."); throw std::runtime_error("Grid size must be positive.");
} }
this->phreeqc = std::make_unique<PhreeqcInit>(database, script); PhreeqcMatrix pqc_mat = PhreeqcMatrix(database, script);
this->phreeqc_mat = pqcScriptToGrid(phreeqc, R); this->transport_names = pqc_mat.getSolutionNames();
this->initial_grid = matToGrid(R, this->phreeqc_mat, grid_def);
const uint32_t solution_count = getSolutionCount(phreeqc, this->initial_grid); Rcpp::Function unique_R("unique");
Rcpp::Function sort_R("sort");
std::vector<std::string> colnames = std::vector<int> unique_ids = Rcpp::as<std::vector<int>>(
Rcpp::as<std::vector<std::string>>(this->initial_grid.names()); unique_R(Rcpp::IntegerVector(grid_def.begin(), grid_def.end())));
this->transport_names = std::vector<std::string>( this->initial_grid = expandGrid(pqc_mat, unique_ids, grid_def);
colnames.begin() + 1,
colnames.begin() + 1 + solution_count); // skip ID
std::map<int, std::string> pqc_raw_dumps; const auto pqc_raw_dumps = replaceRawKeywordIDs(pqc_mat.getDumpStringsPQI());
pqc_raw_dumps = replaceRawKeywordIDs(phreeqc->raw_dumps()); for (const auto &id : unique_ids) {
this->pqc_ids.push_back(id);
this->pqc_ids =
Rcpp::as<std::vector<int>>(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));
} }
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<std::string> colnames =
// Rcpp::as<std::vector<std::string>>(this->initial_grid.names());
// this->transport_names = std::vector<std::string>(
// colnames.begin() + 1,
// colnames.begin() + 1 + solution_count); // skip ID
// std::map<int, std::string> pqc_raw_dumps;
// pqc_raw_dumps = replaceRawKeywordIDs(phreeqc->raw_dumps());
// this->pqc_ids =
// Rcpp::as<std::vector<int>>(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 } // namespace poet

69
src/Init/InitialList.cpp
View File

@ -1,5 +1,6 @@
#include "InitialList.hpp" #include "InitialList.hpp"
#include "DataStructures/NamedVector.hpp" #include "DataStructures/NamedVector.hpp"
#include "PhreeqcMatrix.hpp"
#include <Rcpp/internal/wrap.h> #include <Rcpp/internal/wrap.h>
#include <Rcpp/iostream/Rstreambuf.h> #include <Rcpp/iostream/Rstreambuf.h>
#include <Rcpp/proxy/ProtectedProxy.h> #include <Rcpp/proxy/ProtectedProxy.h>
@ -10,8 +11,8 @@
namespace poet { namespace poet {
void InitialList::initializeFromList(const Rcpp::List &setup) { void InitialList::initializeFromList(const Rcpp::List &setup) {
prepareGrid(setup[grid_key]); PhreeqcMatrix phreeqc = prepareGrid(setup[grid_key]);
initDiffusion(setup[diffusion_key]); initDiffusion(setup[diffusion_key], phreeqc);
initChemistry(setup[chemistry_key]); initChemistry(setup[chemistry_key]);
} }
@ -53,26 +54,26 @@ void InitialList::importList(const Rcpp::List &setup, bool minimal) {
} }
this->database = this->database =
Rcpp::as<std::string>(setup[static_cast<int>(ExportList::CHEM_DATABASE)]); Rcpp::as<std::string>(setup[static_cast<int>(ExportList::CHEM_DATABASE)]);
this->pqc_script = Rcpp::as<std::string>(
setup[static_cast<int>(ExportList::CHEM_PQC_SCRIPT)]);
this->field_header = Rcpp::as<std::vector<std::string>>( this->field_header = Rcpp::as<std::vector<std::string>>(
setup[static_cast<int>(ExportList::CHEM_FIELD_HEADER)]); setup[static_cast<int>(ExportList::CHEM_FIELD_HEADER)]);
this->pqc_scripts = Rcpp::as<std::vector<std::string>>(
setup[static_cast<int>(ExportList::CHEM_PQC_SCRIPTS)]);
this->pqc_ids = Rcpp::as<std::vector<int>>( this->pqc_ids = Rcpp::as<std::vector<int>>(
setup[static_cast<int>(ExportList::CHEM_PQC_IDS)]); setup[static_cast<int>(ExportList::CHEM_PQC_IDS)]);
this->pqc_solutions = Rcpp::as<std::vector<std::string>>( // this->pqc_solutions = Rcpp::as<std::vector<std::string>>(
setup[static_cast<int>(ExportList::CHEM_PQC_SOLUTIONS)]); // setup[static_cast<int>(ExportList::CHEM_PQC_SOLUTIONS)]);
this->pqc_solution_primaries = Rcpp::as<std::vector<std::string>>( // this->pqc_solution_primaries = Rcpp::as<std::vector<std::string>>(
setup[static_cast<int>(ExportList::CHEM_PQC_SOLUTION_PRIMARY)]); // setup[static_cast<int>(ExportList::CHEM_PQC_SOLUTION_PRIMARY)]);
this->pqc_exchanger = // this->pqc_exchanger =
Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_EXCHANGER)]); // Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_EXCHANGER)]);
this->pqc_kinetics = // this->pqc_kinetics =
Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_KINETICS)]); // Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_KINETICS)]);
this->pqc_equilibrium = // this->pqc_equilibrium =
Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_EQUILIBRIUM)]); // Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_EQUILIBRIUM)]);
this->pqc_surface_comps = // this->pqc_surface_comps =
Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_SURFACE_COMPS)]); // Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_SURFACE_COMPS)]);
this->pqc_surface_charges = // this->pqc_surface_charges =
Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_SURFACE_CHARGES)]); // Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_SURFACE_CHARGES)]);
this->dht_species = NamedVector<uint32_t>( this->dht_species = NamedVector<uint32_t>(
setup[static_cast<int>(ExportList::CHEM_DHT_SPECIES)]); setup[static_cast<int>(ExportList::CHEM_DHT_SPECIES)]);
@ -108,25 +109,25 @@ Rcpp::List InitialList::exportList() {
out[static_cast<int>(ExportList::DIFFU_ALPHA_Y)] = this->alpha_y; out[static_cast<int>(ExportList::DIFFU_ALPHA_Y)] = this->alpha_y;
out[static_cast<int>(ExportList::CHEM_DATABASE)] = Rcpp::wrap(this->database); out[static_cast<int>(ExportList::CHEM_DATABASE)] = Rcpp::wrap(this->database);
out[static_cast<int>(ExportList::CHEM_PQC_SCRIPT)] =
Rcpp::wrap(this->pqc_script);
out[static_cast<int>(ExportList::CHEM_FIELD_HEADER)] = out[static_cast<int>(ExportList::CHEM_FIELD_HEADER)] =
Rcpp::wrap(this->field_header); Rcpp::wrap(this->field_header);
out[static_cast<int>(ExportList::CHEM_PQC_SCRIPTS)] =
Rcpp::wrap(this->pqc_scripts);
out[static_cast<int>(ExportList::CHEM_PQC_IDS)] = Rcpp::wrap(this->pqc_ids); out[static_cast<int>(ExportList::CHEM_PQC_IDS)] = Rcpp::wrap(this->pqc_ids);
out[static_cast<int>(ExportList::CHEM_PQC_SOLUTIONS)] = // out[static_cast<int>(ExportList::CHEM_PQC_SOLUTIONS)] =
Rcpp::wrap(this->pqc_solutions); // Rcpp::wrap(this->pqc_solutions);
out[static_cast<int>(ExportList::CHEM_PQC_SOLUTION_PRIMARY)] = // out[static_cast<int>(ExportList::CHEM_PQC_SOLUTION_PRIMARY)] =
Rcpp::wrap(this->pqc_solution_primaries); // Rcpp::wrap(this->pqc_solution_primaries);
out[static_cast<int>(ExportList::CHEM_PQC_EXCHANGER)] = // out[static_cast<int>(ExportList::CHEM_PQC_EXCHANGER)] =
Rcpp::wrap(this->pqc_exchanger); // Rcpp::wrap(this->pqc_exchanger);
out[static_cast<int>(ExportList::CHEM_PQC_KINETICS)] = // out[static_cast<int>(ExportList::CHEM_PQC_KINETICS)] =
Rcpp::wrap(this->pqc_kinetics); // Rcpp::wrap(this->pqc_kinetics);
out[static_cast<int>(ExportList::CHEM_PQC_EQUILIBRIUM)] = // out[static_cast<int>(ExportList::CHEM_PQC_EQUILIBRIUM)] =
Rcpp::wrap(this->pqc_equilibrium); // Rcpp::wrap(this->pqc_equilibrium);
out[static_cast<int>(ExportList::CHEM_PQC_SURFACE_COMPS)] = // out[static_cast<int>(ExportList::CHEM_PQC_SURFACE_COMPS)] =
Rcpp::wrap(this->pqc_surface_comps); // Rcpp::wrap(this->pqc_surface_comps);
out[static_cast<int>(ExportList::CHEM_PQC_SURFACE_CHARGES)] = // out[static_cast<int>(ExportList::CHEM_PQC_SURFACE_CHARGES)] =
Rcpp::wrap(this->pqc_surface_charges); // Rcpp::wrap(this->pqc_surface_charges);
out[static_cast<int>(ExportList::CHEM_DHT_SPECIES)] = this->dht_species; out[static_cast<int>(ExportList::CHEM_DHT_SPECIES)] = this->dht_species;
out[static_cast<int>(ExportList::CHEM_INTERP_SPECIES)] = out[static_cast<int>(ExportList::CHEM_INTERP_SPECIES)] =
Rcpp::wrap(this->interp_species); Rcpp::wrap(this->interp_species);

54
src/Init/InitialList.hpp
View File

@ -2,7 +2,7 @@
#include "Base/RInsidePOET.hpp" #include "Base/RInsidePOET.hpp"
#include "DataStructures/NamedVector.hpp" #include "DataStructures/NamedVector.hpp"
#include "POETInit.hpp" #include "PhreeqcMatrix.hpp"
#include <Rcpp/vector/instantiation.h> #include <Rcpp/vector/instantiation.h>
#include <set> #include <set>
#include <tug/Boundary.hpp> #include <tug/Boundary.hpp>
@ -13,7 +13,6 @@
#include <cstddef> #include <cstddef>
#include <cstdint> #include <cstdint>
#include <PhreeqcInit.hpp>
#include <map> #include <map>
#include <memory> #include <memory>
#include <string> #include <string>
@ -53,22 +52,23 @@ private:
DIFFU_ALPHA_X, DIFFU_ALPHA_X,
DIFFU_ALPHA_Y, DIFFU_ALPHA_Y,
CHEM_DATABASE, CHEM_DATABASE,
CHEM_FIELD_HEADER, CHEM_PQC_SCRIPT,
CHEM_PQC_SCRIPTS,
CHEM_PQC_IDS, CHEM_PQC_IDS,
CHEM_PQC_SOLUTIONS, CHEM_FIELD_HEADER,
CHEM_PQC_SOLUTION_PRIMARY, // CHEM_PQC_SCRIPTS,
CHEM_PQC_EXCHANGER, // CHEM_PQC_SOLUTIONS,
CHEM_PQC_KINETICS, // CHEM_PQC_SOLUTION_PRIMARY,
CHEM_PQC_EQUILIBRIUM, // CHEM_PQC_EXCHANGER,
CHEM_PQC_SURFACE_COMPS, // CHEM_PQC_KINETICS,
CHEM_PQC_SURFACE_CHARGES, // CHEM_PQC_EQUILIBRIUM,
// CHEM_PQC_SURFACE_COMPS,
// CHEM_PQC_SURFACE_CHARGES,
CHEM_DHT_SPECIES, CHEM_DHT_SPECIES,
CHEM_INTERP_SPECIES, CHEM_INTERP_SPECIES,
CHEM_HOOKS, CHEM_HOOKS,
AI_SURROGATE_INPUT_SCRIPT, AI_SURROGATE_INPUT_SCRIPT,
ENUM_SIZE // Hack: Last element of the enum to show enum size ENUM_SIZE // Hack: Last element of the enum to show enum size
}; };
// Grid members // Grid members
static constexpr const char *grid_key = "Grid"; static constexpr const char *grid_key = "Grid";
@ -97,9 +97,9 @@ private:
return GridMembersString[static_cast<std::size_t>(member)]; return GridMembersString[static_cast<std::size_t>(member)];
} }
std::unique_ptr<PhreeqcInit> phreeqc; // std::unique_ptr<PhreeqcMatrix> pqc_mat;
void prepareGrid(const Rcpp::List &grid_input); PhreeqcMatrix prepareGrid(const Rcpp::List &grid_input);
std::uint8_t dim{0}; std::uint8_t dim{0};
@ -114,9 +114,6 @@ private:
Rcpp::List initial_grid; Rcpp::List initial_grid;
// No export
Rcpp::NumericMatrix phreeqc_mat;
public: public:
struct DiffusionInit { struct DiffusionInit {
using BoundaryMap = using BoundaryMap =
@ -169,10 +166,12 @@ private:
return DiffusionMembersString[static_cast<std::size_t>(member)]; return DiffusionMembersString[static_cast<std::size_t>(member)];
} }
void initDiffusion(const Rcpp::List &diffusion_input); void initDiffusion(const Rcpp::List &diffusion_input,
const PhreeqcMatrix &phreeqc);
std::pair<Rcpp::List, Rcpp::List> std::pair<Rcpp::List, Rcpp::List>
resolveBoundaries(const Rcpp::List &boundaries_list, 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 boundaries;
Rcpp::List inner_boundaries; Rcpp::List inner_boundaries;
@ -190,17 +189,10 @@ private:
std::vector<std::string> field_header; std::vector<std::string> field_header;
std::string database; std::string database;
std::vector<std::string> pqc_scripts; std::string pqc_script;
// std::vector<std::string> pqc_scripts;
std::vector<int> pqc_ids; std::vector<int> pqc_ids;
std::vector<std::string> pqc_solutions;
std::vector<std::string> 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<std::uint32_t> dht_species; NamedVector<std::uint32_t> dht_species;
NamedVector<std::uint32_t> interp_species; NamedVector<std::uint32_t> interp_species;
@ -227,10 +219,10 @@ public:
// std::vector<std::string> field_header; // std::vector<std::string> field_header;
std::string database; std::string database;
// std::vector<std::string> pqc_scripts; std::string pqc_script;
// std::vector<int> pqc_ids; std::vector<int> pqc_ids;
std::map<int, POETConfig> pqc_config; // std::map<int, std::string> pqc_input;
// std::vector<std::string> pqc_sol_order; // std::vector<std::string> pqc_sol_order;