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

Add inner_boundaries handling in DiffusionModule and InitialList

Browse Source
This commit is contained in:
Max Luebke 2024-04-05 08:31:32 +00:00
parent 35049e1a77
commit 5cc161e25e
5 changed files with 127 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
ext/tug

@ -1 +1 @@
Subproject commit b104fdcf5238a8a2948b4e835b2d99b64758bb2a
Subproject commit 449647010ab9cdf9e405139f360424a2b21ab3ab

110
src/Init/DiffusionInit.cpp
View File

@ -54,6 +54,7 @@ static std::vector<TugType> colMajToRowMaj(const Rcpp::NumericVector &vec,
static std::vector<std::string> extend_transport_names(
std::unique_ptr<IPhreeqcPOET> &phreeqc, const Rcpp::List &boundaries_list,
const Rcpp::List &inner_boundaries,
const std::vector<std::string> &old_trans_names, Rcpp::List &initial_grid) {
std::vector<std::string> transport_names = old_trans_names;
@ -83,6 +84,13 @@ static std::vector<std::string> extend_transport_names(
}
}
if (inner_boundaries.size() > 0) {
const Rcpp::NumericVector values = inner_boundaries["sol_id"];
for (std::size_t i = 0; i < values.size(); i++) {
constant_pqc_ids.insert(values[i]);
}
}
if (!constant_pqc_ids.empty()) {
for (const auto &pqc_id : constant_pqc_ids) {
const auto solution_names = phreeqc->getSolutionNames(pqc_id);
@ -111,14 +119,17 @@ static Rcpp::List extend_initial_grid(const Rcpp::List &initial_grid,
return extend_grid_R(initial_grid, Rcpp::wrap(names_to_add), old_size);
}
Rcpp::List InitialList::resolveBoundaries(const Rcpp::List &boundaries_list) {
std::pair<Rcpp::List, Rcpp::List>
InitialList::resolveBoundaries(const Rcpp::List &boundaries_list,
const Rcpp::List &inner_boundaries) {
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();
this->transport_names =
extend_transport_names(this->phreeqc, boundaries_list,
extend_transport_names(this->phreeqc, boundaries_list, inner_boundaries,
this->transport_names, this->initial_grid);
const std::size_t new_transport_size = this->transport_names.size();
@ -128,6 +139,10 @@ Rcpp::List InitialList::resolveBoundaries(const Rcpp::List &boundaries_list) {
this->initial_grid, this->transport_names, old_transport_size);
}
const Rcpp::NumericVector &inner_row_vec = inner_boundaries["row"];
const Rcpp::NumericVector &inner_col_vec = inner_boundaries["col"];
const Rcpp::NumericVector &inner_pqc_id_vec = inner_boundaries["sol_id"];
for (const auto &species : this->transport_names) {
Rcpp::List spec_list;
@ -173,9 +188,34 @@ Rcpp::List InitialList::resolveBoundaries(const Rcpp::List &boundaries_list) {
}
bound_list[species] = spec_list;
if (inner_boundaries.size() > 0) {
std::vector<std::uint32_t> rows;
std::vector<std::uint32_t> cols;
std::vector<TugType> c_value;
if (inner_row_vec.size() != inner_col_vec.size() ||
inner_row_vec.size() != inner_pqc_id_vec.size()) {
throw std::runtime_error(
"Inner boundary vectors are not the same length");
}
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), inner_pqc_id_vec[i])));
}
inner_bound[species] =
Rcpp::List::create(Rcpp::Named("row") = Rcpp::wrap(rows),
Rcpp::Named("col") = Rcpp::wrap(cols),
Rcpp::Named("value") = Rcpp::wrap(c_value));
}
}
return bound_list;
return std::make_pair(bound_list, inner_bound);
}
static inline SEXP_TYPE get_datatype(const SEXP &input) {
@ -230,26 +270,36 @@ 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)];
Rcpp::List boundaries;
Rcpp::List inner_boundaries;
if (diffusion_input.containsElementNamed(
DIFFU_MEMBER_STR(DiffusionMembers::BOUNDARIES))) {
boundaries =
diffusion_input[DIFFU_MEMBER_STR(DiffusionMembers::BOUNDARIES)];
}
if (diffusion_input.containsElementNamed(
DIFFU_MEMBER_STR(DiffusionMembers::INNER_BOUNDARIES))) {
inner_boundaries =
diffusion_input[DIFFU_MEMBER_STR(DiffusionMembers::INNER_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);
const auto resolved_boundaries =
resolveBoundaries(boundaries, inner_boundaries);
this->boundaries = resolved_boundaries.first;
this->inner_boundaries = resolved_boundaries.second;
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.parseEval("print(alpha_x)");
// R.parseEval("print(alpha_y)");
}
InitialList::DiffusionInit::BoundaryMap
@ -287,6 +337,39 @@ RcppListToBoundaryMap(const std::vector<std::string> &trans_names,
return map;
}
static InitialList::DiffusionInit::InnerBoundaryMap
RcppListToInnerBoundaryMap(const std::vector<std::string> &trans_names,
const Rcpp::List &inner_bound_list,
std::uint32_t n_cols, std::uint32_t n_rows) {
InitialList::DiffusionInit::InnerBoundaryMap map;
if (inner_bound_list.size() == 0) {
return map;
}
for (const auto &name : trans_names) {
const Rcpp::List &conc_list = inner_bound_list[name];
std::map<std::pair<std::uint32_t, std::uint32_t>, TugType> inner_bc;
const Rcpp::NumericVector &row = conc_list["row"];
const Rcpp::NumericVector &col = conc_list["col"];
const Rcpp::NumericVector &value = conc_list["value"];
if (row.size() != col.size() || row.size() != value.size()) {
throw std::runtime_error(
"Inner boundary vectors are not the same length");
}
for (std::size_t i = 0; i < row.size(); i++) {
inner_bc[std::make_pair(row[i], col[i])] = value[i];
}
map[name] = inner_bc;
}
return map;
}
InitialList::DiffusionInit InitialList::getDiffusionInit() const {
DiffusionInit diff_init;
@ -304,6 +387,9 @@ InitialList::DiffusionInit InitialList::getDiffusionInit() const {
diff_init.boundaries = RcppListToBoundaryMap(
this->transport_names, this->boundaries, this->n_cols, this->n_rows);
diff_init.inner_boundaries =
RcppListToInnerBoundaryMap(this->transport_names, this->inner_boundaries,
this->n_cols, this->n_rows);
diff_init.alpha_x = Field(this->alpha_x);
diff_init.alpha_y = Field(this->alpha_y);

4
src/Init/InitialList.cpp
View File

@ -41,6 +41,8 @@ void InitialList::importList(const Rcpp::List &setup) {
setup[static_cast<int>(ExportList::DIFFU_TRANSPORT)]);
this->boundaries =
Rcpp::List(setup[static_cast<int>(ExportList::DIFFU_BOUNDARIES)]);
this->inner_boundaries =
Rcpp::List(setup[static_cast<int>(ExportList::DIFFU_INNER_BOUNDARIES)]);
this->alpha_x =
Rcpp::List(setup[static_cast<int>(ExportList::DIFFU_ALPHA_X)]);
this->alpha_y =
@ -81,6 +83,8 @@ Rcpp::List InitialList::exportList() {
out[static_cast<int>(ExportList::DIFFU_TRANSPORT)] =
Rcpp::wrap(this->transport_names);
out[static_cast<int>(ExportList::DIFFU_BOUNDARIES)] = this->boundaries;
out[static_cast<int>(ExportList::DIFFU_INNER_BOUNDARIES)] =
this->inner_boundaries;
out[static_cast<int>(ExportList::DIFFU_ALPHA_X)] = this->alpha_x;
out[static_cast<int>(ExportList::DIFFU_ALPHA_Y)] = this->alpha_y;

23
src/Init/InitialList.hpp
View File

@ -16,6 +16,7 @@
#include <IPhreeqcPOET.hpp>
#include <map>
#include <string>
#include <utility>
#include <vector>
#include <DataStructures/Field.hpp>
@ -47,6 +48,7 @@ private:
GRID_INITIAL,
DIFFU_TRANSPORT,
DIFFU_BOUNDARIES,
DIFFU_INNER_BOUNDARIES,
DIFFU_ALPHA_X,
DIFFU_ALPHA_Y,
CHEM_DATABASE,
@ -111,6 +113,10 @@ public:
using BoundaryMap =
std::map<std::string, std::vector<tug::BoundaryElement<TugType>>>;
using InnerBoundaryMap =
std::map<std::string,
std::map<std::pair<std::uint32_t, std::uint32_t>, TugType>>;
uint8_t dim;
std::uint32_t n_cols;
std::uint32_t n_rows;
@ -123,6 +129,7 @@ public:
std::vector<std::string> transport_names;
BoundaryMap boundaries;
InnerBoundaryMap inner_boundaries;
Field alpha_x;
Field alpha_y;
@ -134,22 +141,32 @@ private:
// Diffusion members
static constexpr const char *diffusion_key = "Diffusion";
enum class DiffusionMembers { BOUNDARIES, ALPHA_X, ALPHA_Y, ENUM_SIZE };
enum class DiffusionMembers {
BOUNDARIES,
INNER_BOUNDARIES,
ALPHA_X,
ALPHA_Y,
ENUM_SIZE
};
static constexpr std::size_t size_DiffusionMembers =
static_cast<std::size_t>(InitialList::DiffusionMembers::ENUM_SIZE);
static constexpr std::array<const char *, size_DiffusionMembers>
DiffusionMembersString = {"boundaries", "alpha_x", "alpha_y"};
DiffusionMembersString = {"boundaries", "inner_boundaries", "alpha_x",
"alpha_y"};
constexpr const char *DIFFU_MEMBER_STR(DiffusionMembers member) const {
return DiffusionMembersString[static_cast<std::size_t>(member)];
}
void initDiffusion(const Rcpp::List &diffusion_input);
Rcpp::List resolveBoundaries(const Rcpp::List &boundaries_list);
std::pair<Rcpp::List, Rcpp::List>
resolveBoundaries(const Rcpp::List &boundaries_list,
const Rcpp::List &inner_boundaries);
Rcpp::List boundaries;
Rcpp::List inner_boundaries;
Rcpp::List alpha_x;
Rcpp::List alpha_y;

4
src/Transport/DiffusionModule.cpp
View File

@ -95,6 +95,10 @@ void DiffusionModule::simulate(double requested_dt) {
boundary.deserialize(this->param_list.boundaries[sol_name]);
if (!this->param_list.inner_boundaries[sol_name].empty()) {
boundary.setInnerBoundaries(this->param_list.inner_boundaries[sol_name]);
}
grid.setAlpha(alpha_x, alpha_y);
grid.setConcentrations(conc);