From 5cc161e25ef960a3cae6c96cff5515750465e882 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Fri, 5 Apr 2024 08:31:32 +0000 Subject: [PATCH] Add inner_boundaries handling in DiffusionModule and InitialList --- ext/tug | 2 +- src/Init/DiffusionInit.cpp | 110 ++++++++++++++++++++++++++---- src/Init/InitialList.cpp | 4 ++ src/Init/InitialList.hpp | 23 ++++++- src/Transport/DiffusionModule.cpp | 4 ++ 5 files changed, 127 insertions(+), 16 deletions(-) diff --git a/ext/tug b/ext/tug index b104fdcf5..449647010 160000 --- a/ext/tug +++ b/ext/tug @@ -1 +1 @@ -Subproject commit b104fdcf5238a8a2948b4e835b2d99b64758bb2a +Subproject commit 449647010ab9cdf9e405139f360424a2b21ab3ab diff --git a/src/Init/DiffusionInit.cpp b/src/Init/DiffusionInit.cpp index 6abe43229..8aa8aa96d 100644 --- a/src/Init/DiffusionInit.cpp +++ b/src/Init/DiffusionInit.cpp @@ -54,6 +54,7 @@ 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, Rcpp::List &initial_grid) { std::vector transport_names = old_trans_names; @@ -83,6 +84,13 @@ static std::vector 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 +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 rows; + std::vector cols; + std::vector 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(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 &trans_names, return map; } +static InitialList::DiffusionInit::InnerBoundaryMap +RcppListToInnerBoundaryMap(const std::vector &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, 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); diff --git a/src/Init/InitialList.cpp b/src/Init/InitialList.cpp index 90e0a1dd3..7f5360c87 100644 --- a/src/Init/InitialList.cpp +++ b/src/Init/InitialList.cpp @@ -41,6 +41,8 @@ void InitialList::importList(const Rcpp::List &setup) { setup[static_cast(ExportList::DIFFU_TRANSPORT)]); this->boundaries = Rcpp::List(setup[static_cast(ExportList::DIFFU_BOUNDARIES)]); + this->inner_boundaries = + Rcpp::List(setup[static_cast(ExportList::DIFFU_INNER_BOUNDARIES)]); this->alpha_x = Rcpp::List(setup[static_cast(ExportList::DIFFU_ALPHA_X)]); this->alpha_y = @@ -81,6 +83,8 @@ Rcpp::List InitialList::exportList() { out[static_cast(ExportList::DIFFU_TRANSPORT)] = Rcpp::wrap(this->transport_names); out[static_cast(ExportList::DIFFU_BOUNDARIES)] = this->boundaries; + out[static_cast(ExportList::DIFFU_INNER_BOUNDARIES)] = + this->inner_boundaries; out[static_cast(ExportList::DIFFU_ALPHA_X)] = this->alpha_x; out[static_cast(ExportList::DIFFU_ALPHA_Y)] = this->alpha_y; diff --git a/src/Init/InitialList.hpp b/src/Init/InitialList.hpp index 1500518c1..3b2576ea8 100644 --- a/src/Init/InitialList.hpp +++ b/src/Init/InitialList.hpp @@ -16,6 +16,7 @@ #include #include #include +#include #include #include @@ -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>>; + using InnerBoundaryMap = + std::map, TugType>>; + uint8_t dim; std::uint32_t n_cols; std::uint32_t n_rows; @@ -123,6 +129,7 @@ public: std::vector 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(InitialList::DiffusionMembers::ENUM_SIZE); static constexpr std::array - 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(member)]; } void initDiffusion(const Rcpp::List &diffusion_input); - Rcpp::List resolveBoundaries(const Rcpp::List &boundaries_list); + std::pair + 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; diff --git a/src/Transport/DiffusionModule.cpp b/src/Transport/DiffusionModule.cpp index d81f0d885..2ea80564a 100644 --- a/src/Transport/DiffusionModule.cpp +++ b/src/Transport/DiffusionModule.cpp @@ -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);