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

Add resolvePqcBound function and refactor parseBoundaries2D function

Browse Source
This commit is contained in:
Max Lübke 2024-03-19 15:37:29 +00:00
parent 6ec98077d7
commit bc923b75e6
3 changed files with 85 additions and 49 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
R_lib/init_r_lib.R
View File

@ -28,3 +28,14 @@ pqc_to_grid <- function(pqc_in, grid) {
return(res_df)
}
resolvePqcBound <- function(pqc_mat, transport_spec, id) {
df <- as.data.frame(pqc_mat, check.names = FALSE)
value <- df[df$ID == id, transport_spec]
if (is.nan(value)) {
value <- 0
}
return(value)
}

121
src/Init/DiffusionInit.cpp
View File

@ -26,49 +26,6 @@ const std::map<std::uint8_t, std::string> tug_side_mapping = {
{tug::BC_SIDE_TOP, "N"},
{tug::BC_SIDE_BOTTOM, "S"}};
static Rcpp::List parseBoundaries2D(const Rcpp::List &boundaries_list,
std::uint32_t n_cols,
std::uint32_t n_rows) {
Rcpp::List out_list;
for (const auto &side : tug_side_mapping) {
Rcpp::NumericVector type =
(side.first == tug::BC_SIDE_RIGHT || side.first == tug::BC_SIDE_LEFT)
? Rcpp::NumericVector(n_rows, tug::BC_TYPE_CLOSED)
: Rcpp::NumericVector(n_cols, tug::BC_TYPE_CLOSED);
Rcpp::NumericVector value =
(side.first == tug::BC_SIDE_RIGHT || side.first == tug::BC_SIDE_LEFT)
? Rcpp::NumericVector(n_rows, 0)
: Rcpp::NumericVector(n_cols, 0);
if (boundaries_list.containsElementNamed(side.second.c_str())) {
const Rcpp::List mapping = boundaries_list[side.second];
const Rcpp::NumericVector cells = mapping["cells"];
const Rcpp::NumericVector values = mapping["sol_id"];
if (cells.size() != values.size()) {
throw std::runtime_error("Boundary [" + side.second +
"] cells and values are not the same length");
}
}
out_list[side.second] = Rcpp::List::create(Rcpp::Named("type") = type,
Rcpp::Named("value") = value);
}
}
static inline SEXP_TYPE get_datatype(const SEXP &input) {
Rcpp::Function check_list("is.list");
if (Rcpp::as<bool>(check_list(input))) {
return SEXP_IS_LIST;
} else {
return SEXP_IS_VEC;
}
}
static std::vector<TugType> colMajToRowMaj(const Rcpp::NumericVector &vec,
std::uint32_t n_cols,
std::uint32_t n_rows) {
@ -90,6 +47,70 @@ static std::vector<TugType> colMajToRowMaj(const Rcpp::NumericVector &vec,
}
}
Rcpp::List InitialList::resolveBoundaries(const Rcpp::List &boundaries_list) {
Rcpp::List bound_list;
Rcpp::Function resolve_R("resolvePqcBound");
for (const auto &species : this->transport_names) {
Rcpp::List spec_list;
for (const auto &side : tug_side_mapping) {
Rcpp::NumericVector c_type =
(side.first == tug::BC_SIDE_RIGHT || side.first == tug::BC_SIDE_LEFT)
? Rcpp::NumericVector(n_rows, tug::BC_TYPE_CLOSED)
: Rcpp::NumericVector(n_cols, tug::BC_TYPE_CLOSED);
Rcpp::NumericVector c_value =
(side.first == tug::BC_SIDE_RIGHT || side.first == tug::BC_SIDE_LEFT)
? Rcpp::NumericVector(n_rows, 0)
: Rcpp::NumericVector(n_cols, 0);
if (boundaries_list.containsElementNamed(side.second.c_str())) {
const 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"];
if (cells.size() != values.size()) {
throw std::runtime_error(
"Boundary [" + side.second +
"] cells and values are not the same length");
}
for (std::size_t i = 0; i < cells.size(); i++) {
if (type_str[i] == "closed") {
c_type[cells[i] - 1] = tug::BC_TYPE_CLOSED;
} else if (type_str[i] == "constant") {
c_type[cells[i] - 1] = tug::BC_TYPE_CONSTANT;
c_value[cells[i] - 1] = Rcpp::as<TugType>(
resolve_R(this->phreeqc_mat, Rcpp::wrap(species), values[i]));
} else {
throw std::runtime_error("Unknown boundary type");
}
}
}
spec_list[side.second] = Rcpp::List::create(
Rcpp::Named("type") = c_type, Rcpp::Named("value") = c_value);
}
bound_list[species] = spec_list;
}
return bound_list;
}
static inline SEXP_TYPE get_datatype(const SEXP &input) {
Rcpp::Function check_list("is.list");
if (Rcpp::as<bool>(check_list(input))) {
return SEXP_IS_LIST;
} else {
return SEXP_IS_VEC;
}
}
static Rcpp::List parseAlphas(const SEXP &input,
const std::vector<std::string> &transport_names,
std::uint32_t n_cols, std::uint32_t n_rows) {
@ -132,23 +153,25 @@ static Rcpp::List parseAlphas(const SEXP &input,
return out_list;
}
void InitialList::initDiffusion(const Rcpp::List &diffusion_input) {
// const Rcpp::List &boundaries =
// diffusion_input[DIFFU_MEMBER_STR(DiffusionMembers::BOUNDARIES)];
const Rcpp::List &boundaries =
diffusion_input[DIFFU_MEMBER_STR(DiffusionMembers::BOUNDARIES)];
const SEXP &alpha_x =
diffusion_input[DIFFU_MEMBER_STR(DiffusionMembers::ALPHA_X)];
const SEXP &alpha_y =
diffusion_input[DIFFU_MEMBER_STR(DiffusionMembers::ALPHA_Y)];
this->boundaries = resolveBoundaries(boundaries);
this->alpha_x =
parseAlphas(alpha_x, this->transport_names, this->n_cols, this->n_rows);
this->alpha_y =
parseAlphas(alpha_y, this->transport_names, this->n_cols, this->n_rows);
R["alpha_x"] = this->alpha_x;
R["alpha_y"] = this->alpha_y;
// R["alpha_x"] = this->alpha_x;
// R["alpha_y"] = this->alpha_y;
R.parseEval("print(alpha_x)");
R.parseEval("print(alpha_y)");
// R.parseEval("print(alpha_x)");
// R.parseEval("print(alpha_y)");
}
} // namespace poet

2
src/Init/InitialList.hpp
View File

@ -93,6 +93,8 @@ private:
}
void initDiffusion(const Rcpp::List &diffusion_input);
Rcpp::List resolveBoundaries(const Rcpp::List &boundaries_list);
Rcpp::List boundaries;
Rcpp::List alpha_x;