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

BREAKING CHANGE: Enable Surface/Exchange using new API of PhreeqcEngine

Browse Source
This commit is contained in:
Max Lübke 2024-04-12 19:33:18 +02:00 committed by Max Lübke
parent d7c23c6b3c
commit 480dd146f4
14 changed files with 189 additions and 110 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists.txt
View File

@ -5,7 +5,7 @@ project(POET
DESCRIPTION "A coupled reactive transport simulator")
# specify the C++ standard
set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD 20)
set(CMAKE_CXX_STANDARD_REQUIRED True)
set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE)

2
bench/dolo/dol.pqi
View File

@ -39,5 +39,5 @@ SOLUTION 4
water 1
temp 25
Mg 0.002
Cl 0.002
Cl 0.004
END

4
bench/dolo/dolo_inner_large.R
View File

@ -16,8 +16,8 @@ bound_size <- 2
diffusion_setup <- list(
inner_boundaries = list(
"row" = c(200, 800, 800),
"col" = c(400, 1400, 1600),
"row" = c(400, 1400, 1600),
"col" = c(200, 800, 800),
"sol_id" = c(3, 4, 4)
),
alpha_x = 1e-6,

2
ext/iphreeqc

@ -1 +1 @@
Subproject commit 4ec8b7006c215ad9fc1310505cd56d03ba4b17dd
Subproject commit 749fdbc2e9478046bf3f270c70e5800637246712

17
src/Chemistry/ChemistryModule.cpp
View File

@ -1,12 +1,11 @@
#include "ChemistryModule.hpp"
#include "IPhreeqcPOET.hpp"
#include "PhreeqcEngine.hpp"
#include "SurrogateModels/DHT_Wrapper.hpp"
#include "SurrogateModels/Interpolation.hpp"
#include <algorithm>
#include <cassert>
#include <cstddef>
#include <cstdint>
#include <memory>
#include <mpi.h>
@ -168,12 +167,16 @@ poet::ChemistryModule::ChemistryModule(
this->n_cells = chem_params.total_grid_cells;
if (!is_master) {
for (std::size_t i = 0; i < chem_params.pqc_ids.size(); i++) {
this->phreeqc_instances[chem_params.pqc_ids[i]] =
std::make_unique<IPhreeqcPOET>(chem_params.database,
chem_params.pqc_scripts[i],
chem_params.pqc_sol_order, wp_size_);
for (const auto &[pqc_id, pqc_config] : chem_params.pqc_config) {
this->phreeqc_instances[pqc_id] =
std::make_unique<PhreeqcEngine>(pqc_config);
}
// for (std::size_t i = 0; i < chem_params.pqc_ids.size(); i++) {
// this->phreeqc_instances[chem_params.pqc_ids[i]] =
// std::make_unique<PhreeqcWrapper>(
// chem_params.database, chem_params.pqc_scripts[i],
// chem_params.pqc_sol_order, chem_params.field_header, wp_size_);
// }
}
}

4
src/Chemistry/ChemistryModule.hpp
View File

@ -11,7 +11,7 @@
#include "SurrogateModels/DHT_Wrapper.hpp"
#include "SurrogateModels/Interpolation.hpp"
#include "IPhreeqcPOET.hpp"
#include "PhreeqcEngine.hpp"
#include <array>
#include <cstdint>
#include <map>
@ -378,7 +378,7 @@ protected:
const InitialList::ChemistryInit params;
std::map<int, std::unique_ptr<IPhreeqcPOET>> phreeqc_instances;
std::map<int, std::unique_ptr<PhreeqcEngine>> phreeqc_instances;
};
} // namespace poet

61
src/Chemistry/WorkerFunctions.cpp
View File

@ -4,6 +4,7 @@
#include "Chemistry/ChemistryDefs.hpp"
#include <algorithm>
#include <cstddef>
#include <cstdint>
#include <iomanip>
@ -280,8 +281,6 @@ void poet::ChemistryModule::WorkerReadDHTDump(
void poet::ChemistryModule::WorkerRunWorkPackage(WorkPackage &work_package,
double dSimTime,
double dTimestep) {
// check if we actually need to start phreeqc
std::map<int, std::vector<std::size_t>> zone_mapping;
for (std::size_t wp_id = 0; wp_id < work_package.size; wp_id++) {
if (work_package.mapping[wp_id] != CHEM_PQC) {
@ -289,70 +288,24 @@ void poet::ChemistryModule::WorkerRunWorkPackage(WorkPackage &work_package,
}
auto curr_input = work_package.input[wp_id];
const auto pqc_id = curr_input[0];
const auto pqc_id = static_cast<int>(curr_input[0]);
auto &phreeqc_instance = this->phreeqc_instances[pqc_id];
work_package.output[wp_id] = work_package.input[wp_id];
curr_input.erase(curr_input.begin());
// remove NaNs from the input
curr_input.erase(std::remove_if(curr_input.begin(), curr_input.end(),
[](double d) { return std::isnan(d); }),
curr_input.end());
phreeqc_instance->queueCell(curr_input);
zone_mapping[pqc_id].push_back(wp_id);
}
phreeqc_instance->runCell(curr_input, dTimestep);
if (zone_mapping.empty()) {
return;
}
for (const auto &[pqc_id, eval_cells] : zone_mapping) {
if (eval_cells.empty()) {
continue;
}
this->phreeqc_instances[pqc_id]->runQueuedCells(dTimestep);
std::vector<std::vector<double>> pqc_out;
this->phreeqc_instances[pqc_id]->dequeueCells(pqc_out);
for (std::size_t i = 0; i < eval_cells.size(); i++) {
std::size_t wp_id = eval_cells[i];
std::size_t output_index = 0;
for (std::size_t j = 1; j < work_package.output[wp_id].size(); j++) {
if (!(std::isnan(work_package.output[wp_id][j]))) {
work_package.output[wp_id][j] = pqc_out[i][output_index++];
}
std::size_t output_index = 0;
for (std::size_t i = 0; i < work_package.output[wp_id].size(); i++) {
if (std::isnan(work_package.output[wp_id][i])) {
work_package.output[wp_id][i] = curr_input[output_index++];
}
}
}
// run the phreeqc instances
// for (const auto &[pqc_id, phreeqc_instance] : this->phreeqc_instances) {
// // if (zone_mapping.find(pqc_id) == zone_mapping.end()) {
// // continue;
// // }
// phreeqc_instance->runQueuedCells(dTimestep);
// // remap the output to the work_package
// std::vector<std::vector<double>> pqc_out;
// phreeqc_instance->dequeueCells(pqc_out);
// std::size_t output_id = 0;
// for (const auto &wp_id : zone_mapping[pqc_id]) {
// std::size_t output_index = 0;
// for (std::size_t i = 1; i < work_package.output[wp_id].size(); i++) {
// if (!(std::isnan(work_package.output[wp_id][i]))) {
// work_package.output[wp_id][i] = pqc_out[output_id][output_index++];
// }
// }
// output_id++;
// }
// }
}
void poet::ChemistryModule::WorkerPerfToMaster(int type,

33
src/Init/ChemistryInit.cpp
View File

@ -1,12 +1,16 @@
#include "InitialList.hpp"
#include <Rcpp.h>
#include <cstddef>
#include <vector>
namespace poet {
void InitialList::initChemistry(const Rcpp::List &chem) {
this->pqc_sol_order = this->transport_names;
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"]);
@ -27,6 +31,10 @@ void InitialList::initChemistry(const Rcpp::List &chem) {
}
}
}
this->field_header =
Rcpp::as<std::vector<std::string>>(this->initial_grid.names());
this->field_header.erase(this->field_header.begin());
}
InitialList::ChemistryInit InitialList::getChemistryInit() const {
@ -34,11 +42,26 @@ InitialList::ChemistryInit InitialList::getChemistryInit() const {
chem_init.total_grid_cells = this->n_cols * this->n_rows;
chem_init.database = database;
chem_init.pqc_scripts = pqc_scripts;
chem_init.pqc_ids = pqc_ids;
// chem_init.field_header = this->field_header;
chem_init.pqc_sol_order = pqc_sol_order;
chem_init.database = database;
// 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};
}
// chem_init.pqc_sol_order = pqc_solutions;
chem_init.dht_species = dht_species;
chem_init.interp_species = interp_species;

6
src/Init/DiffusionInit.cpp
View File

@ -7,8 +7,8 @@
#include <set>
#include "DataStructures/Field.hpp"
#include "IPhreeqcPOET.hpp"
#include "InitialList.hpp"
#include "PhreeqcInit.hpp"
#include <Rcpp/Function.h>
#include <Rcpp/proxy/ProtectedProxy.h>
@ -54,7 +54,7 @@ static std::vector<TugType> colMajToRowMaj(const Rcpp::NumericVector &vec,
}
static std::vector<std::string>
extend_transport_names(std::unique_ptr<IPhreeqcPOET> &phreeqc,
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) {
@ -332,7 +332,7 @@ RcppListToBoundaryMap(const std::vector<std::string> &trans_names,
}
for (std::size_t i = 0; i < type.size(); i++) {
if (type[i] == tug::BC_TYPE_CONSTANT) {
if (static_cast<tug::BC_TYPE>(type[i]) == tug::BC_TYPE_CONSTANT) {
bc.setBoundaryElementConstant(static_cast<tug::BC_SIDE>(tug_index), i,
value[i]);
}

18
src/Init/GridInit.cpp
View File

@ -1,6 +1,7 @@
#include "InitialList.hpp"
#include <IPhreeqcPOET.hpp>
#include "PhreeqcInit.hpp"
#include <RInside.h>
#include <Rcpp/Function.h>
#include <Rcpp/vector/instantiation.h>
@ -15,8 +16,8 @@
namespace poet {
static Rcpp::NumericMatrix
pqcScriptToGrid(std::unique_ptr<IPhreeqcPOET> &phreeqc, RInside &R) {
IPhreeqcPOET::PhreeqcMat phreeqc_mat = phreeqc->getPhreeqcMat();
pqcScriptToGrid(std::unique_ptr<PhreeqcInit> &phreeqc, RInside &R) {
PhreeqcInit::PhreeqcMat phreeqc_mat = phreeqc->getPhreeqcMat();
// add "id" to the front of the column names
@ -59,9 +60,9 @@ replaceRawKeywordIDs(std::map<int, std::string> raws) {
return raws;
}
static inline uint32_t getSolutionCount(std::unique_ptr<IPhreeqcPOET> &phreeqc,
static inline uint32_t getSolutionCount(std::unique_ptr<PhreeqcInit> &phreeqc,
const Rcpp::List &initial_grid) {
IPhreeqcPOET::ModulesArray mod_array;
PhreeqcInit::ModulesArray mod_array;
Rcpp::Function unique_R("unique");
std::vector<int> row_ids =
@ -155,7 +156,7 @@ void InitialList::prepareGrid(const Rcpp::List &grid_input) {
throw std::runtime_error("Grid size must be positive.");
}
this->phreeqc = std::make_unique<IPhreeqcPOET>(database, script);
this->phreeqc = std::make_unique<PhreeqcInit>(database, script);
this->phreeqc_mat = pqcScriptToGrid(phreeqc, R);
this->initial_grid = matToGrid(R, this->phreeqc_mat, grid_def);
@ -178,6 +179,11 @@ void InitialList::prepareGrid(const Rcpp::List &grid_input) {
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));
}
}

36
src/Init/InitialList.cpp
View File

@ -53,12 +53,26 @@ void InitialList::importList(const Rcpp::List &setup, bool minimal) {
}
this->database =
Rcpp::as<std::string>(setup[static_cast<int>(ExportList::CHEM_DATABASE)]);
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_sol_order = Rcpp::as<std::vector<std::string>>(
setup[static_cast<int>(ExportList::CHEM_PQC_SOL_ORDER)]);
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)]);
@ -92,11 +106,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_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_SOL_ORDER)] =
Rcpp::wrap(this->pqc_sol_order);
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);

37
src/Init/InitialList.hpp
View File

@ -2,8 +2,8 @@
#include "Base/RInsidePOET.hpp"
#include "DataStructures/NamedVector.hpp"
#include "POETInit.hpp"
#include <Rcpp/vector/instantiation.h>
#include <memory>
#include <set>
#include <tug/Boundary.hpp>
@ -13,8 +13,9 @@
#include <cstddef>
#include <cstdint>
#include <IPhreeqcPOET.hpp>
#include <PhreeqcInit.hpp>
#include <map>
#include <memory>
#include <string>
#include <utility>
#include <vector>
@ -52,9 +53,16 @@ private:
DIFFU_ALPHA_X,
DIFFU_ALPHA_Y,
CHEM_DATABASE,
CHEM_FIELD_HEADER,
CHEM_PQC_SCRIPTS,
CHEM_PQC_IDS,
CHEM_PQC_SOL_ORDER,
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,
@ -88,7 +96,7 @@ private:
return GridMembersString[static_cast<std::size_t>(member)];
}
std::unique_ptr<IPhreeqcPOET> phreeqc;
std::unique_ptr<PhreeqcInit> phreeqc;
void prepareGrid(const Rcpp::List &grid_input);
@ -178,11 +186,19 @@ private:
void initChemistry(const Rcpp::List &chem_input);
std::vector<std::string> field_header;
std::string database;
std::vector<std::string> pqc_scripts;
std::vector<int> pqc_ids;
std::vector<std::string> pqc_sol_order;
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;
@ -204,10 +220,15 @@ public:
struct ChemistryInit {
uint32_t total_grid_cells;
// std::vector<std::string> field_header;
std::string database;
std::vector<std::string> pqc_scripts;
std::vector<int> pqc_ids;
std::vector<std::string> pqc_sol_order;
// std::vector<std::string> pqc_scripts;
// std::vector<int> pqc_ids;
std::map<int, POETConfig> pqc_config;
// std::vector<std::string> pqc_sol_order;
NamedVector<std::uint32_t> dht_species;
NamedVector<std::uint32_t> interp_species;

46
src/poet.cpp
View File

@ -65,15 +65,16 @@ static void init_global_functions(RInside &R) {
// HACK: this is a step back as the order and also the count of fields is
// predefined, but it will change in the future
void writeFieldsToR(RInside &R, const Field &trans, const Field &chem) {
// static inline void writeFieldsToR(RInside &R, const Field &trans,
// const Field &chem) {
Rcpp::DataFrame t_field(trans.asSEXP());
R["TMP"] = t_field;
R.parseEval("mysetup$state_T <- TMP");
// Rcpp::DataFrame t_field(trans.asSEXP());
// R["TMP"] = t_field;
// R.parseEval("mysetup$state_T <- TMP");
R["TMP"] = chem.asSEXP();
R.parseEval("mysetup$state_C <- TMP");
}
// R["TMP"] = chem.asSEXP();
// R.parseEval("mysetup$state_C <- TMP");
// }
enum ParseRet { PARSER_OK, PARSER_ERROR, PARSER_HELP };
@ -222,7 +223,26 @@ ParseRet parseInitValues(char **argv, RuntimeParameters &params) {
return ParseRet::PARSER_OK;
}
static Rcpp::List RunMasterLoop(const RuntimeParameters &params,
// HACK: this is a step back as the order and also the count of fields is
// predefined, but it will change in the future
void call_master_iter_end(RInside &R, const Field &trans, const Field &chem) {
R["TMP"] = Rcpp::wrap(trans.AsVector());
R["TMP_PROPS"] = Rcpp::wrap(trans.GetProps());
R.parseEval(std::string("state_T <- setNames(data.frame(matrix(TMP, nrow=" +
std::to_string(trans.GetRequestedVecSize()) +
")), TMP_PROPS)"));
R["TMP"] = Rcpp::wrap(chem.AsVector());
R["TMP_PROPS"] = Rcpp::wrap(chem.GetProps());
R.parseEval(std::string("state_C <- setNames(data.frame(matrix(TMP, nrow=" +
std::to_string(chem.GetRequestedVecSize()) +
")), TMP_PROPS)"));
R["setup"] = *global_rt_setup;
R.parseEval("setup <- master_iteration_end(setup, state_T, state_C)");
*global_rt_setup = R["setup"];
}
static Rcpp::List RunMasterLoop(RInsidePOET &R, const RuntimeParameters &params,
DiffusionModule &diffusion,
ChemistryModule &chem) {
@ -263,13 +283,7 @@ static Rcpp::List RunMasterLoop(const RuntimeParameters &params,
// state_C after every iteration if the cmdline option
// --ignore-results is not given (and thus the R variable
// store_result is TRUE)
{
const auto &trans_field = diffusion.getField().asSEXP();
const auto &chem_field = chem.getField().asSEXP();
*global_rt_setup = master_iteration_end_R.value()(
*global_rt_setup, trans_field, chem_field);
}
call_master_iter_end(R, diffusion.getField(), chem.getField());
diffusion.getField().update(chem.getField());
@ -400,7 +414,7 @@ int main(int argc, char *argv[]) {
chemistry.masterSetField(init_list.getInitialGrid());
Rcpp::List profiling = RunMasterLoop(run_params, diffusion, chemistry);
Rcpp::List profiling = RunMasterLoop(R, run_params, diffusion, chemistry);
MSG("finished simulation loop");

31
util/data_evaluation/jlFun_Eval.jl Normal file
View File

@ -0,0 +1,31 @@
using RData, DataFrames, Plots
# Load all .rds files in given directory
function load_data(directory)
files = readdir(directory)
# data is dictionary with iteration number as key
data = Dict{Int,Any}()
for file in files
if (endswith(file, ".rds") && startswith(file, "iter"))
# extract iteration number (iter_XXX.rds)
iter = parse(Int, split(split(file, "_")[2], ".")[1])
data[iter] = load(joinpath(directory, file))
end
end
return data
end
function spec_to_mat(in_df::DataFrame, spec::Symbol, cols::Int)
specvec = in_df[!, spec]
specmat = transpose(reshape(specvec, cols, :))
return specmat
end
function plot_field(data::AbstractMatrix, log::Bool=false)
if log
heatmap(1:size(data, 2), 1:size(data, 1), log10.(data), c=:viridis)
else
heatmap(1:size(data, 2), 1:size(data, 1), data, c=:viridis)
end
end