mirror of
https://git.gfz-potsdam.de/naaice/poet.git
synced 2025-12-15 12:28:22 +01:00
Update subproject commit in ext/iphreeqc
This commit is contained in:
Max Lübke
parent
6316c43094
commit
45ea77ae0f
@ -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
|
||||
|
||||
@ -1 +1 @@
|
||||
Subproject commit e6e5e0d5156c093241a53e6ce074ef346d64ae26
|
||||
Subproject commit 61214a01ad4cf99527f657e6a41c68282e4886c4
|
||||
@ -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<PhreeqcEngine>(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<PhreeqcEngine>(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]] =
|
||||
|
||||
@ -2,16 +2,14 @@
|
||||
|
||||
#include <Rcpp.h>
|
||||
#include <cstddef>
|
||||
#include <fstream>
|
||||
#include <string>
|
||||
#include <utility>
|
||||
#include <vector>
|
||||
|
||||
namespace poet {
|
||||
|
||||
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")) {
|
||||
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.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<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};
|
||||
}
|
||||
// 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;
|
||||
|
||||
|
||||
@ -1,4 +1,3 @@
|
||||
#include <algorithm>
|
||||
#include <tug/Boundary.hpp>
|
||||
// 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 <Rcpp/Function.h>
|
||||
#include <Rcpp/proxy/ProtectedProxy.h>
|
||||
@ -53,98 +52,102 @@ static std::vector<TugType> colMajToRowMaj(const Rcpp::NumericVector &vec,
|
||||
}
|
||||
}
|
||||
|
||||
static std::vector<std::string>
|
||||
extend_transport_names(std::unique_ptr<PhreeqcInit> &phreeqc,
|
||||
const Rcpp::List &boundaries_list,
|
||||
const Rcpp::List &inner_boundaries,
|
||||
const std::vector<std::string> &old_trans_names) {
|
||||
// static std::vector<std::string>
|
||||
// extend_transport_names(std::unique_ptr<PhreeqcInit> &phreeqc,
|
||||
// const Rcpp::List &boundaries_list,
|
||||
// const Rcpp::List &inner_boundaries,
|
||||
// const std::vector<std::string> &old_trans_names) {
|
||||
|
||||
std::vector<std::string> transport_names = old_trans_names;
|
||||
std::set<int> constant_pqc_ids;
|
||||
// std::vector<std::string> transport_names = old_trans_names;
|
||||
// std::set<int> 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<std::size_t>(values[i]));
|
||||
}
|
||||
}
|
||||
}
|
||||
// for (auto i = 0; i < cells.size(); i++) {
|
||||
// if (type_str[i] == "constant") {
|
||||
// constant_pqc_ids.insert(static_cast<std::size_t>(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<std::size_t>(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<std::size_t>(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<std::string> &transport_names) {
|
||||
// std::vector<std::string> 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<std::string> &transport_names) {
|
||||
// // std::vector<std::string> 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<Rcpp::List, Rcpp::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 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<TugType>(
|
||||
resolve_R(this->phreeqc_mat, Rcpp::wrap(species_name), values[i]));
|
||||
c_value[c_id] =
|
||||
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 {
|
||||
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<TugType>(resolve_R(
|
||||
this->phreeqc_mat, Rcpp::wrap(species_name), inner_pqc_id_vec[i])));
|
||||
c_value.push_back(static_cast<TugType>(
|
||||
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] =
|
||||
@ -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;
|
||||
|
||||
|
||||
@ -1,11 +1,13 @@
|
||||
#include "InitialList.hpp"
|
||||
|
||||
#include "PhreeqcInit.hpp"
|
||||
#include "PhreeqcMatrix.hpp"
|
||||
|
||||
#include <RInside.h>
|
||||
#include <Rcpp/Function.h>
|
||||
#include <Rcpp/vector/Matrix.h>
|
||||
#include <Rcpp/vector/instantiation.h>
|
||||
#include <cstdint>
|
||||
#include <fstream>
|
||||
#include <map>
|
||||
#include <memory>
|
||||
#include <regex>
|
||||
@ -15,37 +17,41 @@
|
||||
|
||||
namespace poet {
|
||||
|
||||
static Rcpp::NumericMatrix
|
||||
pqcScriptToGrid(std::unique_ptr<PhreeqcInit> &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<int, std::string>
|
||||
replaceRawKeywordIDs(std::map<int, std::string> raws) {
|
||||
@ -60,24 +66,25 @@ replaceRawKeywordIDs(std::map<int, std::string> raws) {
|
||||
return raws;
|
||||
}
|
||||
|
||||
static inline uint32_t getSolutionCount(std::unique_ptr<PhreeqcInit> &phreeqc,
|
||||
const Rcpp::List &initial_grid) {
|
||||
PhreeqcInit::ModulesArray mod_array;
|
||||
Rcpp::Function unique_R("unique");
|
||||
// static inline uint32_t getSolutionCount(std::unique_ptr<PhreeqcInit>
|
||||
// &phreeqc,
|
||||
// const Rcpp::List &initial_grid) {
|
||||
// PhreeqcInit::ModulesArray mod_array;
|
||||
// Rcpp::Function unique_R("unique");
|
||||
|
||||
std::vector<int> row_ids =
|
||||
Rcpp::as<std::vector<int>>(unique_R(initial_grid["ID"]));
|
||||
// std::vector<int> row_ids =
|
||||
// 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");
|
||||
// sizes_vec = Rcpp::as<std::vector<std::uint32_t>>(
|
||||
// modify_sizes(sizes_vec, phreeqc_mat, initial_grid));
|
||||
// // Rcpp::Function modify_sizes("modify_module_sizes");
|
||||
// // sizes_vec = Rcpp::as<std::vector<std::uint32_t>>(
|
||||
// // 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<int> 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<PhreeqcInit>(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<std::string> colnames =
|
||||
Rcpp::as<std::vector<std::string>>(this->initial_grid.names());
|
||||
std::vector<int> unique_ids = Rcpp::as<std::vector<int>>(
|
||||
unique_R(Rcpp::IntegerVector(grid_def.begin(), grid_def.end())));
|
||||
|
||||
this->transport_names = std::vector<std::string>(
|
||||
colnames.begin() + 1,
|
||||
colnames.begin() + 1 + solution_count); // skip ID
|
||||
this->initial_grid = expandGrid(pqc_mat, unique_ids, grid_def);
|
||||
|
||||
std::map<int, std::string> pqc_raw_dumps;
|
||||
const auto pqc_raw_dumps = replaceRawKeywordIDs(pqc_mat.getDumpStringsPQI());
|
||||
|
||||
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));
|
||||
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<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
|
||||
@ -1,5 +1,6 @@
|
||||
#include "InitialList.hpp"
|
||||
#include "DataStructures/NamedVector.hpp"
|
||||
#include "PhreeqcMatrix.hpp"
|
||||
#include <Rcpp/internal/wrap.h>
|
||||
#include <Rcpp/iostream/Rstreambuf.h>
|
||||
#include <Rcpp/proxy/ProtectedProxy.h>
|
||||
@ -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<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>>(
|
||||
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>>(
|
||||
setup[static_cast<int>(ExportList::CHEM_PQC_IDS)]);
|
||||
this->pqc_solutions = Rcpp::as<std::vector<std::string>>(
|
||||
setup[static_cast<int>(ExportList::CHEM_PQC_SOLUTIONS)]);
|
||||
this->pqc_solution_primaries = Rcpp::as<std::vector<std::string>>(
|
||||
setup[static_cast<int>(ExportList::CHEM_PQC_SOLUTION_PRIMARY)]);
|
||||
this->pqc_exchanger =
|
||||
Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_EXCHANGER)]);
|
||||
this->pqc_kinetics =
|
||||
Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_KINETICS)]);
|
||||
this->pqc_equilibrium =
|
||||
Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_EQUILIBRIUM)]);
|
||||
this->pqc_surface_comps =
|
||||
Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_SURFACE_COMPS)]);
|
||||
this->pqc_surface_charges =
|
||||
Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_SURFACE_CHARGES)]);
|
||||
// this->pqc_solutions = Rcpp::as<std::vector<std::string>>(
|
||||
// setup[static_cast<int>(ExportList::CHEM_PQC_SOLUTIONS)]);
|
||||
// this->pqc_solution_primaries = Rcpp::as<std::vector<std::string>>(
|
||||
// setup[static_cast<int>(ExportList::CHEM_PQC_SOLUTION_PRIMARY)]);
|
||||
// this->pqc_exchanger =
|
||||
// Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_EXCHANGER)]);
|
||||
// this->pqc_kinetics =
|
||||
// Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_KINETICS)]);
|
||||
// this->pqc_equilibrium =
|
||||
// Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_EQUILIBRIUM)]);
|
||||
// this->pqc_surface_comps =
|
||||
// Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_SURFACE_COMPS)]);
|
||||
// this->pqc_surface_charges =
|
||||
// Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_SURFACE_CHARGES)]);
|
||||
|
||||
this->dht_species = NamedVector<uint32_t>(
|
||||
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::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)] =
|
||||
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_SOLUTIONS)] =
|
||||
Rcpp::wrap(this->pqc_solutions);
|
||||
out[static_cast<int>(ExportList::CHEM_PQC_SOLUTION_PRIMARY)] =
|
||||
Rcpp::wrap(this->pqc_solution_primaries);
|
||||
out[static_cast<int>(ExportList::CHEM_PQC_EXCHANGER)] =
|
||||
Rcpp::wrap(this->pqc_exchanger);
|
||||
out[static_cast<int>(ExportList::CHEM_PQC_KINETICS)] =
|
||||
Rcpp::wrap(this->pqc_kinetics);
|
||||
out[static_cast<int>(ExportList::CHEM_PQC_EQUILIBRIUM)] =
|
||||
Rcpp::wrap(this->pqc_equilibrium);
|
||||
out[static_cast<int>(ExportList::CHEM_PQC_SURFACE_COMPS)] =
|
||||
Rcpp::wrap(this->pqc_surface_comps);
|
||||
out[static_cast<int>(ExportList::CHEM_PQC_SURFACE_CHARGES)] =
|
||||
Rcpp::wrap(this->pqc_surface_charges);
|
||||
// out[static_cast<int>(ExportList::CHEM_PQC_SOLUTIONS)] =
|
||||
// Rcpp::wrap(this->pqc_solutions);
|
||||
// out[static_cast<int>(ExportList::CHEM_PQC_SOLUTION_PRIMARY)] =
|
||||
// Rcpp::wrap(this->pqc_solution_primaries);
|
||||
// out[static_cast<int>(ExportList::CHEM_PQC_EXCHANGER)] =
|
||||
// Rcpp::wrap(this->pqc_exchanger);
|
||||
// out[static_cast<int>(ExportList::CHEM_PQC_KINETICS)] =
|
||||
// Rcpp::wrap(this->pqc_kinetics);
|
||||
// out[static_cast<int>(ExportList::CHEM_PQC_EQUILIBRIUM)] =
|
||||
// Rcpp::wrap(this->pqc_equilibrium);
|
||||
// out[static_cast<int>(ExportList::CHEM_PQC_SURFACE_COMPS)] =
|
||||
// Rcpp::wrap(this->pqc_surface_comps);
|
||||
// out[static_cast<int>(ExportList::CHEM_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_INTERP_SPECIES)] =
|
||||
Rcpp::wrap(this->interp_species);
|
||||
|
||||
@ -2,7 +2,7 @@
|
||||
|
||||
#include "Base/RInsidePOET.hpp"
|
||||
#include "DataStructures/NamedVector.hpp"
|
||||
#include "POETInit.hpp"
|
||||
#include "PhreeqcMatrix.hpp"
|
||||
#include <Rcpp/vector/instantiation.h>
|
||||
#include <set>
|
||||
#include <tug/Boundary.hpp>
|
||||
@ -13,7 +13,6 @@
|
||||
#include <cstddef>
|
||||
#include <cstdint>
|
||||
|
||||
#include <PhreeqcInit.hpp>
|
||||
#include <map>
|
||||
#include <memory>
|
||||
#include <string>
|
||||
@ -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<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};
|
||||
|
||||
@ -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<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>
|
||||
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<std::string> field_header;
|
||||
|
||||
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<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> interp_species;
|
||||
@ -227,10 +219,10 @@ public:
|
||||
// std::vector<std::string> field_header;
|
||||
|
||||
std::string database;
|
||||
// std::vector<std::string> pqc_scripts;
|
||||
// std::vector<int> pqc_ids;
|
||||
std::string pqc_script;
|
||||
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;
|
||||
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user