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

Add Basic heterogeneous diffusion functionality

Browse Source
This commit is contained in:
Max Lübke 2024-03-27 15:35:53 +00:00
parent a04344e435
commit 0bfc95bb52
16 changed files with 889 additions and 1083 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

148
R_lib/kin_r_library.R
View File

@ -28,7 +28,7 @@ master_init <- function(setup) {
## Setup the directory where we will store the results
verb <- FALSE
if (local_rank == 0) {
# if (local_rank == 0) {
verb <- TRUE ## verbosity loading MUFITS results
if (!dir.exists(fileout)) {
dir.create(fileout)
@ -41,25 +41,26 @@ master_init <- function(setup) {
} else {
msgm("store_result is ", store_result)
}
} else {
}
# } else {
# }
setup$iter <- 1
setup$maxiter <- setup$iterations
setup$timesteps <- setup$timesteps
setup$maxiter <- length(setup$timesteps)
setup$iterations <- setup$maxiter
setup$simulation_time <- 0
if (is.null(setup[["store_result"]])) {
setup$store_result <- TRUE
}
if (setup$store_result) {
if (is.null(setup[["out_save"]])) {
setup$out_save <- seq(1, setup$iterations)
}
}
return(setup)
}
@ -67,31 +68,33 @@ master_init <- function(setup) {
## calculated time step if store_result is TRUE and increments the
## iteration counter
master_iteration_end <- function(setup) {
iter <- setup$iter
# iter <- setup$iter
# print(iter)
## max digits for iterations
dgts <- as.integer(ceiling(log10(setup$iterations + 1)))
## string format to use in sprintf
iter <- 1
dgts <- as.integer(ceiling(log10(1)))
## string format to use in sprintf
fmt <- paste0("%0", dgts, "d")
## Write on disk state_T and state_C after every iteration
## comprised in setup$out_save
if (setup$store_result) {
if (iter %in% setup$out_save) {
nameout <- paste0(fileout, "/iter_", sprintf(fmt=fmt, iter), ".rds")
info <- list(
tr_req_dt = as.integer(setup$req_dt)
## tr_allow_dt = setup$allowed_dt,
## tr_inniter = as.integer(setup$inniter)
)
saveRDS(list(
T = setup$state_T, C = setup$state_C,
simtime = as.integer(setup$simtime),
tr_info = info
), file = nameout)
msgm("results stored in <", nameout, ">")
}
}
msgm("done iteration", iter, "/", setup$maxiter)
# if (setup$store_result) {
# if (iter %in% setup$out_save) {
nameout <- paste0(fileout, "/iter_", sprintf(fmt = fmt, iter), ".rds")
info <- list(
tr_req_dt = as.integer(1)
## tr_allow_dt = setup$allowed_dt,
## tr_inniter = as.integer(setup$inniter)
)
saveRDS(list(
T = setup$state_T, C = setup$state_C,
simtime = as.integer(0),
tr_info = info
), file = nameout)
msgm("results stored in <", nameout, ">")
# }
# }
msgm("done iteration", iter, "/", 1)
setup$iter <- setup$iter + 1
return(setup)
}
@ -104,41 +107,41 @@ slave_chemistry <- function(setup, data) {
immobile <- setup$immobile
kin <- setup$kin
ann <- setup$ann
iter <- setup$iter
timesteps <- setup$timesteps
dt <- timesteps[iter]
state_T <- data ## not the global field, but the work-package
## treat special H+/pH, e-/pe cases
state_T <- RedModRphree::Act2pH(state_T)
## reduction of the problem
if (setup$reduce) {
reduced <- ReduceStateOmit(state_T, omit = setup$ann)
} else {
reduced <- state_T
}
## form the PHREEQC input script for the current work package
inplist <- SplitMultiKin(
data = reduced, procs = 1, base = base, first = first,
ann = ann, prop = prop, minerals = immobile, kin = kin, dt = dt
)
## if (local_rank==1 & iter==1)
## RPhreeWriteInp("FirstInp", inplist)
tmpC <- RunPQC(inplist, procs = 1, second = TRUE)
## recompose after the reduction
if (setup$reduce) {
state_C <- RecomposeState(tmpC, reduced)
} else {
state_C <- tmpC
}
## the next line is needed since we don't need all columns of
## PHREEQC output
return(state_C[, prop])
@ -147,40 +150,40 @@ slave_chemistry <- function(setup, data) {
## This function, called by master
master_chemistry <- function(setup, data) {
state_T <- setup$state_T
msgm(" chemistry iteration", setup$iter)
## treat special H+/pH, e-/pe cases
state_T <- RedModRphree::Act2pH(state_T)
## reduction of the problem
if (setup$reduce) {
reduced <- ReduceStateOmit(state_T, omit = setup$ann)
} else {
reduced <- state_T
}
## inject data from workers
res_C <- data
rownames(res_C) <- NULL
## print(res_C)
if (nrow(res_C) > nrow(reduced)) {
res_C <- res_C[seq(2, nrow(res_C), by = 2), ]
}
## recompose after the reduction
if (setup$reduce) {
state_C <- RecomposeState(res_C, reduced)
} else {
state_C <- res_C
}
setup$state_C <- state_C
setup$reduced <- reduced
return(setup)
}
@ -188,7 +191,7 @@ master_chemistry <- function(setup, data) {
## Adapted version for "reduction"
ReduceStateOmit <- function(data, omit = NULL, sign = 6) {
require(mgcv)
rem <- colnames(data)
if (is.list(omit)) {
indomi <- match(names(omit), colnames(data))
@ -196,13 +199,13 @@ ReduceStateOmit <- function(data, omit = NULL, sign = 6) {
} else {
datao <- data
}
datao <- signif(datao, sign)
red <- mgcv::uniquecombs(datao)
inds <- attr(red, "index")
now <- ncol(red)
## reattach the omitted column(s)
## FIXME: control if more than one ann is present
if (is.list(omit)) {
@ -224,11 +227,11 @@ ReduceStateOmit <- function(data, omit = NULL, sign = 6) {
## Attach the name of the calling function to the message displayed on
## R's stdout
msgm <- function(...) {
if (local_rank == 0) {
fname <- as.list(sys.call(-1))[[1]]
prefix <- paste0("R: ", fname, " ::")
cat(paste(prefix, ..., "\n"))
}
# if (local_rank == 0) {
fname <- as.list(sys.call(-1))[[1]]
prefix <- paste0("R: ", fname, " ::")
cat(paste(prefix, ..., "\n"))
# }
invisible()
}
@ -236,18 +239,17 @@ msgm <- function(...) {
## Function called by master R process to store on disk all relevant
## parameters for the simulation
StoreSetup <- function(setup) {
to_store <- vector(mode = "list", length = 4)
## names(to_store) <- c("Sim", "Flow", "Transport", "Chemistry", "DHT")
names(to_store) <- c("Sim", "Transport", "DHT", "Cmdline")
## read the setup R file, which is sourced in kin.cpp
tmpbuff <- file(filesim, "r")
setupfile <- readLines(tmpbuff)
close.connection(tmpbuff)
to_store$Sim <- setupfile
## to_store$Flow <- list(
## snapshots = setup$snapshots,
## gridfile = setup$gridfile,
@ -282,22 +284,22 @@ StoreSetup <- function(setup) {
to_store$DHT <- FALSE
}
if (dht_enabled) {
to_store$DHT <- list(
enabled = dht_enabled,
log = dht_log
#signif = dht_final_signif,
#proptype = dht_final_proptype
)
} else {
to_store$DHT <- FALSE
}
if (dht_enabled) {
to_store$DHT <- list(
enabled = dht_enabled,
log = dht_log
# signif = dht_final_signif,
# proptype = dht_final_proptype
)
} else {
to_store$DHT <- FALSE
}
saveRDS(to_store, file = paste0(fileout, "/setup.rds"))
msgm("initialization stored in ", paste0(fileout, "/setup.rds"))
saveRDS(to_store, file = paste0(fileout, "/setup.rds"))
msgm("initialization stored in ", paste0(fileout, "/setup.rds"))
}
GetWorkPackageSizesVector <- function(n_packages, package_size, len) {
ids <- rep(1:n_packages, times=package_size, each = 1)[1:len]
ids <- rep(1:n_packages, times = package_size, each = 1)[1:len]
return(as.integer(table(ids)))
}

4
bench/barite/het/runtime.R Normal file
View File

@ -0,0 +1,4 @@
list(
timesteps = c(1),
store_result = TRUE
)

2
ext/tug

@ -1 +1 @@
Subproject commit 25855da6b2930559b542bbadb16299932332d6a3
Subproject commit 8c0687a6cd4a10a79c7a554083a35eda11cc55f0

274
src/Base/Grid.cpp
View File

@ -1,274 +0,0 @@
/*
** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
** Potsdam)
**
** Copyright (C) 2018-2023 Marco De Lucia, Max Luebke (GFZ Potsdam)
**
** POET is free software; you can redistribute it and/or modify it under the
** terms of the GNU General Public License as published by the Free Software
** Foundation; either version 2 of the License, or (at your option) any later
** version.
**
** POET is distributed in the hope that it will be useful, but WITHOUT ANY
** WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
** A PARTICULAR PURPOSE. See the GNU General Public License for more details.
**
** You should have received a copy of the GNU General Public License along with
** this program; if not, write to the Free Software Foundation, Inc., 51
** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
#include "Grid.hpp"
#include "SimParams.hpp"
#include <RInside.h>
#include <Rcpp.h>
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <new>
#include <numeric>
#include <stdexcept>
#include <string>
#include <vector>
using namespace poet;
using namespace Rcpp;
enum { INIT_SCRATCH, INIT_PHREEQC, INIT_RDS };
static inline int8_t get_type_id(std::string type) {
if (type == "scratch") {
return INIT_SCRATCH;
}
if (type == "phreeqc") {
return INIT_PHREEQC;
}
if (type == "rds") {
return INIT_RDS;
}
return -1;
}
void poet::Grid::InitModuleFromParams(const poet::GridParams &grid_args) {
assert(grid_args.n_cells.size() == grid_args.s_cells.size());
this->SetGridDimension(grid_args.n_cells.size());
std::copy(grid_args.n_cells.begin(), grid_args.n_cells.end(),
this->n_cells.begin());
std::copy(grid_args.s_cells.begin(), grid_args.s_cells.end(),
this->grid_size.begin());
auto init_type = get_type_id(grid_args.type);
assert(init_type >= 0);
switch (init_type) {
case INIT_SCRATCH:
this->SetPropNames(grid_args.props);
this->InitGridFromScratch(grid_args.init_df);
break;
case INIT_RDS:
this->SetPropNames(grid_args.props);
this->InitGridFromRDS();
break;
case INIT_PHREEQC:
this->InitGridFromPhreeqc();
}
}
void poet::Grid::SetGridDimension(uint8_t dim) {
assert(dim < MAX_DIM + 1);
this->dim = dim;
}
void poet::Grid::SetGridCellCount(uint32_t n_x, uint32_t n_y, uint32_t n_z) {
assert(this->dim < MAX_DIM + 1);
this->n_cells[0] = n_x;
this->n_cells[1] = n_y;
// this->n_cells[2] = n_z;
}
void poet::Grid::SetGridSize(double s_x, double s_y, double s_z) {
assert(this->dim < MAX_DIM + 1);
this->grid_size[0] = s_x;
this->grid_size[1] = s_y;
// this->grid_size[2] = s_z;
}
void poet::Grid::SetPropNames(const std::vector<std::string> &prop_names) {
this->prop_names = prop_names;
}
void poet::Grid::PushbackModuleFlow(const std::string &input,
const std::string &output) {
FlowInputOutputInfo element = {input, output};
this->flow_vec.push_back(element);
}
void poet::Grid::PreModuleFieldCopy(uint32_t tick) {
FlowInputOutputInfo curr_element = this->flow_vec.at(tick);
const std::string input_module_name = curr_element.input_field;
StateMemory *out_state = this->GetStatePointer(curr_element.output_field);
std::vector<std::string> &mod_props = out_state->props;
std::vector<double> &mod_field = out_state->mem;
// copy output of another module as input for another module
for (uint32_t i = 0; i < mod_props.size(); i++) {
try {
std::vector<double> t_prop_vec =
this->GetSpeciesByName(mod_props[i], input_module_name);
std::copy(t_prop_vec.begin(), t_prop_vec.end(),
mod_field.begin() + (i * this->GetTotalCellCount()));
} catch (...) {
// TODO: there might be something wrong ...
continue;
}
}
}
Grid::Grid() {
this->n_cells.fill(0);
this->grid_size.fill(0);
};
Grid::~Grid() {
for (auto &[key, val] : this->state_register) {
delete val;
}
}
void poet::Grid::InitGridFromScratch(const Rcpp::DataFrame &init_cell) {
const uint32_t vec_size = this->GetTotalCellCount();
StateMemory *curr_state =
this->RegisterState(poet::GRID_MODULE_NAME, this->prop_names);
curr_state->props = this->prop_names;
std::vector<double> &curr_field = curr_state->mem;
for (auto const &prop : this->prop_names) {
// size of the vector shall be 1
double prop_val = ((Rcpp::NumericVector)init_cell[prop.c_str()])[0];
std::vector<double> prop_vec(vec_size, prop_val);
curr_field.insert(curr_field.end(), prop_vec.begin(), prop_vec.end());
}
this->grid_init = curr_state;
}
void poet::Grid::InitGridFromRDS() {
// TODO
}
void poet::Grid::InitGridFromPhreeqc() {
// TODO
}
poet::StateMemory *Grid::RegisterState(std::string module_name,
std::vector<std::string> props) {
poet::StateMemory *mem = new poet::StateMemory;
mem->props = props;
const auto [it, success] = this->state_register.insert({module_name, mem});
if (!success) {
delete mem;
throw std::bad_alloc();
}
return mem;
}
poet::StateMemory *Grid::GetStatePointer(std::string module_name) {
const auto it = this->state_register.find(module_name);
if (it == this->state_register.end()) {
throw std::range_error(
std::string("Requested module " + module_name + " was not found"));
}
return it->second;
}
auto Grid::GetGridSize(uint8_t direction) const -> uint32_t {
return this->grid_size.at(direction);
}
auto Grid::GetGridDimension() const -> uint8_t { return this->dim; }
auto Grid::GetTotalCellCount() const -> uint32_t {
uint32_t sum = 1;
for (auto const &dim : this->n_cells) {
sum *= (dim != 0 ? dim : 1);
}
return sum;
}
auto Grid::GetGridCellsCount(uint8_t direction) const -> uint32_t {
return this->n_cells.at(direction);
}
auto poet::Grid::GetInitialGrid() const -> StateMemory * {
return this->grid_init;
}
auto Grid::GetSpeciesCount() const -> uint32_t {
return this->prop_names.size();
}
auto Grid::GetPropNames() const -> std::vector<std::string> {
return this->prop_names;
}
auto poet::Grid::GetSpeciesByName(std::string name,
std::string module_name) const
-> std::vector<double> {
auto const it = this->state_register.find(module_name);
if (it == this->state_register.end()) {
throw std::range_error(
std::string("Requested module " + module_name + " was not found"));
}
poet::StateMemory const *module_memory = it->second;
auto const &prop_vector = module_memory->props;
auto const find_res = std::find(prop_vector.begin(), prop_vector.end(), name);
if (find_res == prop_vector.end()) {
throw std::range_error(std::string(
"Species " + name + " was not found for module " + module_name));
}
uint32_t prop_index = find_res - prop_vector.begin();
uint32_t begin_vec = prop_index * this->GetTotalCellCount(),
end_vec = ((prop_index + 1) * this->GetTotalCellCount());
return std::vector<double>(module_memory->mem.begin() + begin_vec,
module_memory->mem.begin() + end_vec);
}
// HACK: Helper function
void poet::Grid::WriteFieldsToR(RInside &R) {
for (auto const &elem : this->state_register) {
std::string field_name = elem.first;
StateMemory *field = elem.second;
const uint32_t n_cells_per_prop = field->mem.size() / field->props.size();
R["TMP"] = Rcpp::wrap(field->mem);
R["TMP_PROPS"] = Rcpp::wrap(field->props);
R.parseEval(std::string(
"mysetup$" + field_name + " <- setNames(data.frame(matrix(TMP, nrow=" +
std::to_string(n_cells_per_prop) + ")), TMP_PROPS)"));
}
}

117
src/Base/Grid.hpp
View File

@ -1,117 +0,0 @@
/*
** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
** Potsdam)
**
** Copyright (C) 2018-2023 Marco De Lucia, Max Luebke (GFZ Potsdam)
**
** POET is free software; you can redistribute it and/or modify it under the
** terms of the GNU General Public License as published by the Free Software
** Foundation; either version 2 of the License, or (at your option) any later
** version.
**
** POET is distributed in the hope that it will be useful, but WITHOUT ANY
** WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
** A PARTICULAR PURPOSE. See the GNU General Public License for more details.
**
** You should have received a copy of the GNU General Public License along with
** this program; if not, write to the Free Software Foundation, Inc., 51
** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
#ifndef GRID_H
#define GRID_H
#include "SimParams.hpp"
#include <RInside.h>
#include <Rcpp.h>
#include <array>
#include <cstddef>
#include <cstdint>
#include <list>
#include <map>
#include <string>
#include <vector>
#define MAX_DIM 2
namespace poet {
enum { GRID_X_DIR, GRID_Y_DIR, GRID_Z_DIR };
using StateMemory = struct s_StateMemory {
std::vector<double> mem;
std::vector<std::string> props;
};
using FlowInputOutputInfo = struct inOut_info {
std::string input_field;
std::string output_field;
};
constexpr const char *GRID_MODULE_NAME = "grid_init";
/**
* @brief Class describing the grid
*
* Providing methods to shuffle and unshuffle grid (for the master) as also to
* import and export a work package (for worker).
*
* @todo find better abstraction
*
*/
class Grid {
public:
Grid();
~Grid();
void InitModuleFromParams(const poet::GridParams &grid_args);
void SetGridDimension(uint8_t dim);
void SetGridCellCount(uint32_t n_x, uint32_t n_y = 0, uint32_t n_z = 0);
void SetGridSize(double s_x, double s_y = 0., double s_z = 0.);
void SetPropNames(const std::vector<std::string> &prop_names);
void PushbackModuleFlow(const std::string &input, const std::string &output);
void PreModuleFieldCopy(uint32_t tick);
void InitGridFromScratch(const Rcpp::DataFrame &init_cell);
void InitGridFromRDS();
void InitGridFromPhreeqc();
auto GetGridDimension() const -> uint8_t;
auto GetTotalCellCount() const -> uint32_t;
auto GetGridCellsCount(uint8_t direction) const -> uint32_t;
auto GetGridSize(uint8_t direction) const -> uint32_t;
auto RegisterState(std::string module_name, std::vector<std::string> props)
-> StateMemory *;
auto GetStatePointer(std::string module_name) -> StateMemory *;
auto GetInitialGrid() const -> StateMemory *;
auto GetSpeciesCount() const -> uint32_t;
auto GetPropNames() const -> std::vector<std::string>;
auto GetSpeciesByName(std::string name,
std::string module_name = poet::GRID_MODULE_NAME) const
-> std::vector<double>;
void WriteFieldsToR(RInside &R);
private:
std::vector<FlowInputOutputInfo> flow_vec;
std::uint8_t dim = 0;
std::array<double, MAX_DIM> grid_size;
std::array<std::uint32_t, MAX_DIM> n_cells;
std::map<std::string, StateMemory *> state_register;
StateMemory *grid_init = std::nullptr_t();
std::vector<std::string> prop_names;
};
} // namespace poet
#endif // GRID_H

163
src/Base/SimParams.cpp
View File

@ -1,10 +1,11 @@
// Time-stamp: "Last modified 2023-08-15 12:12:36 delucia"
/*
** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
** Potsdam)
**
** Copyright (C) 2018-2023 Marco De Lucia, Max Luebke (GFZ Potsdam)
** Copyright (C) 2018-2022 Marco De Lucia, Max Luebke (GFZ Potsdam)
**
** Copyright (C) 2018-2022 Marco De Lucia (GFZ Potsdam), Max Luebke (University
** of Potsdam)
**
** POET is free software; you can redistribute it and/or modify it under the
** terms of the GNU General Public License as published by the Free Software
@ -22,104 +23,33 @@
#include "SimParams.hpp"
#include "../Chemistry/enums.hpp"
#include "Base/Macros.hpp"
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <RInside.h>
#include <Rcpp.h>
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <string>
#include <string_view>
#include <vector>
#include "Base/RInsidePOET.hpp"
#include "argh.hpp" // Argument handler https://github.com/adishavit/argh
using namespace poet;
using namespace std;
using namespace Rcpp;
poet::GridParams::s_GridParams(RInside &R) {
auto tmp_n_cells =
Rcpp::as<std::vector<uint32_t>>(R.parseEval("mysetup$grid$n_cells"));
assert(tmp_n_cells.size() < 3);
RuntimeParameters::RuntimeParameters(RInsidePOET &R, char *argv[],
int world_rank_)
: world_rank(world_rank_) {
this->dim = tmp_n_cells.size();
std::copy(tmp_n_cells.begin(), tmp_n_cells.end(), this->n_cells.begin());
auto tmp_s_cells =
Rcpp::as<std::vector<double>>(R.parseEval("mysetup$grid$s_cells"));
assert(tmp_s_cells.size() == this->dim);
std::copy(tmp_s_cells.begin(), tmp_s_cells.end(), this->s_cells.begin());
this->total_n =
(dim == 1 ? this->n_cells[0] : this->n_cells[0] * this->n_cells[1]);
this->type = Rcpp::as<std::string>(R.parseEval("mysetup$grid$type"));
}
poet::DiffusionParams::s_DiffusionParams(RInside &R) {
this->initial_t =
Rcpp::as<Rcpp::DataFrame>(R.parseEval("mysetup$diffusion$init"));
this->alpha =
Rcpp::as<Rcpp::NumericVector>(R.parseEval("mysetup$diffusion$alpha"));
if (Rcpp::as<bool>(
R.parseEval("'vecinj_inner' %in% names(mysetup$diffusion)"))) {
this->vecinj_inner =
Rcpp::as<Rcpp::List>(R.parseEval("mysetup$diffusion$vecinj_inner"));
}
this->vecinj =
Rcpp::as<Rcpp::DataFrame>(R.parseEval("mysetup$diffusion$vecinj"));
this->vecinj_index =
Rcpp::as<Rcpp::DataFrame>(R.parseEval("mysetup$diffusion$vecinj_index"));
}
void poet::ChemistryParams::initFromR(RInsidePOET &R) {
this->database_path =
Rcpp::as<std::string>(R.parseEval("mysetup$chemistry$database"));
this->input_script =
Rcpp::as<std::string>(R.parseEval("mysetup$chemistry$input_script"));
if (R.checkIfExists("dht_species", "mysetup$chemistry")) {
this->dht_signifs = Rcpp::as<NamedVector<std::uint32_t>>(
R.parseEval(("mysetup$chemistry$dht_species")));
}
if (R.checkIfExists("pht_species", "mysetup$chemistry")) {
this->pht_signifs = Rcpp::as<NamedVector<std::uint32_t>>(
R.parseEval(("mysetup$chemistry$pht_species")));
}
this->hooks.dht_fill =
RHookFunction<bool>(R, "mysetup$chemistry$hooks$dht_fill");
this->hooks.dht_fuzz =
RHookFunction<std::vector<double>>(R, "mysetup$chemistry$hooks$dht_fuzz");
this->hooks.interp_pre = RHookFunction<std::vector<std::size_t>>(
R, "mysetup$chemistry$hooks$interp_pre_func");
this->hooks.interp_post =
RHookFunction<bool>(R, "mysetup$chemistry$hooks$interp_post_func");
}
SimParams::SimParams(int world_rank_, int world_size_) {
this->simparams.world_rank = world_rank_;
this->simparams.world_size = world_size_;
}
int SimParams::parseFromCmdl(char *argv[], RInsidePOET &R) {
// initialize argh object
argh::parser cmdl(argv);
// if user asked for help
if (cmdl[{"help", "h"}]) {
if (simparams.world_rank == 0) {
if (this->world_rank == 0) {
MSG("Todo");
MSG("See README.md for further information.");
}
return poet::PARSER_HELP;
}
// if positional arguments are missing
@ -242,11 +172,74 @@ int SimParams::parseFromCmdl(char *argv[], RInsidePOET &R) {
this->chem_params.initFromR(R);
return poet::PARSER_OK;
};
// poet::GridParams::s_GridParams(RInside &R) {
// auto tmp_n_cells =
// Rcpp::as<std::vector<uint32_t>>(R.parseEval("mysetup$grid$n_cells"));
// assert(tmp_n_cells.size() < 3);
// this->dim = tmp_n_cells.size();
// std::copy(tmp_n_cells.begin(), tmp_n_cells.end(), this->n_cells.begin());
// auto tmp_s_cells =
// Rcpp::as<std::vector<double>>(R.parseEval("mysetup$grid$s_cells"));
// assert(tmp_s_cells.size() == this->dim);
// std::copy(tmp_s_cells.begin(), tmp_s_cells.end(), this->s_cells.begin());
// this->total_n =
// (dim == 1 ? this->n_cells[0] : this->n_cells[0] * this->n_cells[1]);
// this->type = Rcpp::as<std::string>(R.parseEval("mysetup$grid$type"));
// }
// poet::DiffusionParams::s_DiffusionParams(RInside &R) {
// this->initial_t =
// Rcpp::as<Rcpp::DataFrame>(R.parseEval("mysetup$diffusion$init"));
// this->alpha =
// Rcpp::as<Rcpp::NumericVector>(R.parseEval("mysetup$diffusion$alpha"));
// if (Rcpp::as<bool>(
// R.parseEval("'vecinj_inner' %in% names(mysetup$diffusion)"))) {
// this->vecinj_inner =
// Rcpp::as<Rcpp::List>(R.parseEval("mysetup$diffusion$vecinj_inner"));
// }
// this->vecinj =
// Rcpp::as<Rcpp::DataFrame>(R.parseEval("mysetup$diffusion$vecinj"));
// this->vecinj_index =
// Rcpp::as<Rcpp::DataFrame>(R.parseEval("mysetup$diffusion$vecinj_index"));
// }
void poet::ChemistryParams::initFromR(RInsidePOET &R) {
// this->database_path =
// Rcpp::as<std::string>(R.parseEval("mysetup$chemistry$database"));
// this->input_script =
// Rcpp::as<std::string>(R.parseEval("mysetup$chemistry$input_script"));
if (R.checkIfExists("dht_species", "mysetup$chemistry")) {
this->dht_signifs = Rcpp::as<NamedVector<std::uint32_t>>(
R.parseEval(("mysetup$chemistry$dht_species")));
}
if (R.checkIfExists("pht_species", "mysetup$chemistry")) {
this->pht_signifs = Rcpp::as<NamedVector<std::uint32_t>>(
R.parseEval(("mysetup$chemistry$pht_species")));
}
this->hooks.dht_fill =
RHookFunction<bool>(R, "mysetup$chemistry$hooks$dht_fill");
this->hooks.dht_fuzz =
RHookFunction<std::vector<double>>(R, "mysetup$chemistry$hooks$dht_fuzz");
this->hooks.interp_pre = RHookFunction<std::vector<std::size_t>>(
R, "mysetup$chemistry$hooks$interp_pre_func");
this->hooks.interp_post =
RHookFunction<bool>(R, "mysetup$chemistry$hooks$interp_post_func");
}
void SimParams::initVectorParams(RInside &R) {}
void RuntimeParameters::initVectorParams(RInside &R) {}
std::list<std::string> SimParams::validateOptions(argh::parser cmdl) {
std::list<std::string> RuntimeParameters::validateOptions(argh::parser cmdl) {
/* store all unknown parameters here */
std::list<std::string> retList;

299
src/Base/SimParams.hpp
View File

@ -4,6 +4,9 @@
**
** Copyright (C) 2018-2022 Marco De Lucia, Max Luebke (GFZ Potsdam)
**
** Copyright (C) 2018-2022 Marco De Lucia (GFZ Potsdam), Max Luebke (University
** of Potsdam)
**
** POET is free software; you can redistribute it and/or modify it under the
** terms of the GNU General Public License as published by the Free Software
** Foundation; either version 2 of the License, or (at your option) any later
@ -18,20 +21,13 @@
** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
#ifndef PARSER_H
#define PARSER_H
#pragma once
#include "../DataStructures/DataStructures.hpp"
#include "Macros.hpp"
#include "DataStructures/NamedVector.hpp"
#include "RInsidePOET.hpp"
#include "argh.hpp" // Argument handler https://github.com/adishavit/argh
#include <cstdint>
#include <optional>
#include <string>
#include <string_view>
#include <unordered_map>
#include <vector>
#include <RInside.h>
@ -51,7 +47,8 @@ enum { PARSER_OK, PARSER_ERROR, PARSER_HELP };
* @brief Defining all simulation parameters
*
*/
typedef struct {
struct RuntimeParameters {
/** Count of processes in MPI_COMM_WORLD */
int world_size;
/** rank of proces in MPI_COMM_WORLD */
@ -79,40 +76,40 @@ typedef struct {
bool print_progressbar;
bool interp_enabled;
} t_simparams;
using GridParams = struct s_GridParams {
std::array<uint32_t, 2> n_cells;
std::array<double, 2> s_cells;
std::uint8_t dim;
std::uint32_t total_n;
std::string type;
Rcpp::DataFrame init_df;
std::string input_script;
std::string database_path;
std::vector<std::string> props;
s_GridParams(RInside &R);
};
using DiffusionParams = struct s_DiffusionParams {
Rcpp::DataFrame initial_t;
// using GridParams = struct s_GridParams {
// std::array<uint32_t, 2> n_cells;
// std::array<double, 2> s_cells;
// std::uint8_t dim;
// std::uint32_t total_n;
Rcpp::NumericVector alpha;
Rcpp::List vecinj_inner;
// std::string type;
Rcpp::DataFrame vecinj;
Rcpp::DataFrame vecinj_index;
// Rcpp::DataFrame init_df;
// std::string input_script;
// std::string database_path;
s_DiffusionParams(RInside &R);
};
// std::vector<std::string> props;
// s_GridParams(RInside &R);
// };
// using DiffusionParams = struct s_DiffusionParams {
// Rcpp::DataFrame initial_t;
// Rcpp::NumericVector alpha;
// Rcpp::List vecinj_inner;
// Rcpp::DataFrame vecinj;
// Rcpp::DataFrame vecinj_index;
// s_DiffusionParams(RInside &R);
// };
struct ChemistryParams {
std::string database_path;
std::string input_script;
// std::string database_path;
// std::string input_script;
bool use_dht;
std::uint64_t dht_size;
@ -137,6 +134,7 @@ struct ChemistryParams {
void initFromR(RInsidePOET &R);
};
} // namespace poet
/**
* @brief Reads information from program arguments and R runtime
*
@ -147,129 +145,136 @@ struct ChemistryParams {
* Stores and distribute current simulation parameters at any time.
*
*/
class SimParams {
public:
/**
* @brief Construct a new SimParams object
*
* With all given parameters a new instance of this class will be created.
*
* @param world_rank Rank of process inside MPI_COMM_WORLD
* @param world_size Size of communicator MPI_COMM_WORLD
*/
SimParams(int world_rank, int world_size);
// class SimParams {
// public:
// /**
// * @brief Construct a new SimParams object
// *
// * With all given parameters a new instance of this class will be created.
// *
// * @param world_rank Rank of process inside MPI_COMM_WORLD
// * @param world_size Size of communicator MPI_COMM_WORLD
// */
// SimParams(int world_rank);
/**
* @brief Parse program arguments
*
* This is done by the argh.h library.
*
* First, the function will check if there is a flag 'help' or 'h'. If this is
* the case a help message is printed and the function will return with
* PARSER_HELP.
*
* Second, if there are not 2 positional arguments an error will be printed to
* stderr and the function returns with PARSER_ERROR.
*
* Then all given program parameters and flags will be read and checked, if
* there are known by validateOptions. A list of all unknown options might be
* returned, printed out and the function will return with PARSER_ERROR.
* Oterhwise the function continuos.
*
* Now all program arguments will be stored inside t_simparams struct, printed
* out and the function returns with PARSER_OK.
*
* Also, all parsed agruments are distributed to the R runtime.
*
* @param argv Argument value of the program
* @param R Instantiated R runtime
* @return int Returns with 0 if no error occured, otherwise value less than 0
* is returned.
*/
int parseFromCmdl(char *argv[], RInsidePOET &R);
// /**
// * @brief Parse program arguments
// *
// * This is done by the argh.h library.
// *
// * First, the function will check if there is a flag 'help' or 'h'. If this
// is
// * the case a help message is printed and the function will return with
// * PARSER_HELP.
// *
// * Second, if there are not 2 positional arguments an error will be printed
// to
// * stderr and the function returns with PARSER_ERROR.
// *
// * Then all given program parameters and flags will be read and checked, if
// * there are known by validateOptions. A list of all unknown options might
// be
// * returned, printed out and the function will return with PARSER_ERROR.
// * Oterhwise the function continuos.
// *
// * Now all program arguments will be stored inside t_simparams struct,
// printed
// * out and the function returns with PARSER_OK.
// *
// * Also, all parsed agruments are distributed to the R runtime.
// *
// * @param argv Argument value of the program
// * @param R Instantiated R runtime
// * @return int Returns with 0 if no error occured, otherwise value less
// than 0
// * is returned.
// */
// int parseFromCmdl(char *argv[], RInsidePOET &R);
/**
* @brief Init std::vector values
*
* This will initialize dht_signif_vector and dht_prop_type_vector internally
* depending on whether vectors are defined by R-Simulation file or not.
* 'init_chemistry' must be run beforehand.
*
* @param R R runtime
*/
void initVectorParams(RInside &R);
// /**
// * @brief Init std::vector values
// *
// * This will initialize dht_signif_vector and dht_prop_type_vector
// internally
// * depending on whether vectors are defined by R-Simulation file or not.
// * 'init_chemistry' must be run beforehand.
// *
// * @param R R runtime
// */
// void initVectorParams(RInside &R);
/**
* @brief Get the numerical params struct
*
* Returns a struct which contains all numerical or boolean simulation
* parameters.
*
* @return t_simparams Parameter struct
*/
auto getNumParams() const { return this->simparams; };
// /**
// * @brief Get the numerical params struct
// *
// * Returns a struct which contains all numerical or boolean simulation
// * parameters.
// *
// * @return t_simparams Parameter struct
// */
// auto getNumParams() const { return this->simparams; };
/**
* @brief Get the DHT_Signif_Vector
*
* Returns a vector indicating which significant values are used for each
* variable of a grid cell.
*
* @return std::vector<int> Vector of integers containing information about
* significant digits
*/
auto getDHTSignifVector() const { return this->dht_signif_vector; };
// /**
// * @brief Get the DHT_Signif_Vector
// *
// * Returns a vector indicating which significant values are used for each
// * variable of a grid cell.
// *
// * @return std::vector<int> Vector of integers containing information about
// * significant digits
// */
// auto getDHTSignifVector() const { return this->dht_signif_vector; };
auto getPHTSignifVector() const { return this->pht_signif_vector; };
// auto getPHTSignifVector() const { return this->pht_signif_vector; };
auto getPHTBucketSize() const { return this->pht_bucket_size; };
auto getPHTMinEntriesNeeded() const { return this->pht_min_entries_needed; };
// auto getPHTBucketSize() const { return this->pht_bucket_size; };
// auto getPHTMinEntriesNeeded() const { return this->pht_min_entries_needed;
// };
/**
* @brief Get the filesim name
*
* Returns a string containing the absolute path to a R file defining the
* simulation.
*
* @return std::string Absolute path to R file
*/
auto getFilesim() const { return this->filesim; };
// /**
// * @brief Get the filesim name
// *
// * Returns a string containing the absolute path to a R file defining the
// * simulation.
// *
// * @return std::string Absolute path to R file
// */
// auto getFilesim() const { return this->filesim; };
/**
* @brief Get the output directory
*
* Returns the name of an absolute path where all output files should be
* stored.
*
* @return std::string Absolute path to output path
*/
auto getOutDir() const { return this->out_dir; };
// /**
// * @brief Get the output directory
// *
// * Returns the name of an absolute path where all output files should be
// * stored.
// *
// * @return std::string Absolute path to output path
// */
// auto getOutDir() const { return this->out_dir; };
const auto &getChemParams() const { return chem_params; }
// const auto &getChemParams() const { return chem_params; }
private:
std::list<std::string> validateOptions(argh::parser cmdl);
// private:
// std::list<std::string> validateOptions(argh::parser cmdl);
const std::set<std::string> flaglist{"ignore-result", "dht", "P", "progress",
"interp"};
const std::set<std::string> paramlist{
"work-package-size", "dht-strategy",
"dht-size", "dht-snaps",
"dht-file", "interp-size",
"interp-min", "interp-bucket-entries"};
// const std::set<std::string> flaglist{"ignore-result", "dht", "P",
// "progress",
// "interp"};
// const std::set<std::string> paramlist{
// "work-package-size", "dht-strategy",
// "dht-size", "dht-snaps",
// "dht-file", "interp-size",
// "interp-min", "interp-bucket-entries"};
t_simparams simparams;
// t_simparams simparams;
std::vector<uint32_t> dht_signif_vector;
std::vector<uint32_t> pht_signif_vector;
// std::vector<uint32_t> dht_signif_vector;
// std::vector<uint32_t> pht_signif_vector;
uint32_t pht_bucket_size;
uint32_t pht_min_entries_needed;
// uint32_t pht_bucket_size;
// uint32_t pht_min_entries_needed;
std::string filesim;
std::string out_dir;
// std::string filesim;
// std::string out_dir;
ChemistryParams chem_params;
};
} // namespace poet
#endif // PARSER_H
// ChemistryParams chem_params;
// };
// } // namespace poet

17
src/CMakeLists.txt
View File

@ -1,14 +1,13 @@
add_library(POETLib OBJECT)
add_library(POETLib)
target_sources(
POETLib
target_sources(POETLib
PRIVATE
Init/InitialList.cpp
Init/GridInit.cpp
Init/DiffusionInit.cpp
Init/ChemistryInit.cpp
PRIVATE
DataStructures/Field.cpp
Transport/DiffusionModule.cpp
)
target_include_directories(POETLib PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}")
@ -16,7 +15,7 @@ target_link_libraries(
POETLib
PUBLIC RRuntime
PUBLIC IPhreeqcPOET
PRIVATE tug
PUBLIC tug
)
# add_library(poetlib
@ -43,7 +42,7 @@ target_link_libraries(
# tug
# )
# target_compile_definitions(poetlib PUBLIC STRICT_R_HEADERS OMPI_SKIP_MPICXX)
target_compile_definitions(POETLib PUBLIC STRICT_R_HEADERS OMPI_SKIP_MPICXX)
# mark_as_advanced(PHREEQCRM_BUILD_MPI PHREEQCRM_DISABLE_OPENMP)
# set(PHREEQCRM_DISABLE_OPENMP ON CACHE BOOL "" FORCE)
@ -65,9 +64,9 @@ file(READ "${PROJECT_SOURCE_DIR}/R_lib/init_r_lib.R" R_INIT_LIB)
configure_file(poet.hpp.in poet.hpp @ONLY)
# add_executable(poet poet.cpp)
# target_link_libraries(poet PRIVATE poetlib MPI::MPI_C RRuntime)
# target_include_directories(poet PRIVATE "${CMAKE_CURRENT_BINARY_DIR}")
add_executable(poet poet.cpp)
target_link_libraries(poet PRIVATE POETLib MPI::MPI_C RRuntime)
target_include_directories(poet PRIVATE "${CMAKE_CURRENT_BINARY_DIR}")
add_executable(poet_init initializer.cpp)
target_link_libraries(poet_init PRIVATE POETLib RRuntime)

6
src/Init/CMakeLists.txt
View File

@ -1,6 +0,0 @@
target_sources(POETLib PRIVATE
InitialList.cpp
GridInit.cpp
DiffusionInit.cpp
ChemistryInit.cpp
)

42
src/Init/DiffusionInit.cpp
View File

@ -5,6 +5,7 @@
#include <map>
#include <set>
#include "DataStructures/Field.hpp"
#include "InitialList.hpp"
#include <Rcpp/proxy/ProtectedProxy.h>
@ -14,6 +15,7 @@
#include <cstdint>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
namespace poet {
@ -175,6 +177,42 @@ void InitialList::initDiffusion(const Rcpp::List &diffusion_input) {
// R.parseEval("print(alpha_y)");
}
InitialList::DiffusionInit::BoundaryMap
RcppListToBoundaryMap(const std::vector<std::string> &trans_names,
const Rcpp::List &bound_list, uint32_t n_cols,
uint32_t n_rows) {
InitialList::DiffusionInit::BoundaryMap map;
for (const auto &name : trans_names) {
const Rcpp::List &conc_list = bound_list[name];
tug::Boundary<TugType> bc(n_rows, n_cols);
for (const auto &[tug_index, r_init_name] : tug_side_mapping) {
const Rcpp::List &side_list = conc_list[tug_index];
const Rcpp::NumericVector &type = side_list["type"];
const Rcpp::NumericVector &value = side_list["value"];
if (type.size() != value.size()) {
throw std::runtime_error(
"Boundary type and value are not the same length");
}
for (std::size_t i = 0; i < type.size(); i++) {
if (type[i] == tug::BC_TYPE_CONSTANT) {
bc.setBoundaryElementConstant(static_cast<tug::BC_SIDE>(tug_index), i,
value[i]);
}
}
}
map[name] = bc.serialize();
}
return map;
}
InitialList::DiffusionInit InitialList::getDiffusionInit() const {
DiffusionInit diff_init;
@ -189,8 +227,8 @@ InitialList::DiffusionInit InitialList::getDiffusionInit() const {
diff_init.constant_cells = this->constant_cells;
diff_init.transport_names = this->transport_names;
diff_init.initial_grid = Field(this->initial_grid);
diff_init.boundaries = Field(this->boundaries);
diff_init.boundaries = RcppListToBoundaryMap(
this->transport_names, this->boundaries, this->n_cols, this->n_rows);
diff_init.alpha_x = Field(this->alpha_x);
diff_init.alpha_y = Field(this->alpha_y);

8
src/Init/InitialList.cpp
View File

@ -38,10 +38,10 @@ 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->alpha_x = Rcpp::as<std::vector<double>>(
setup[static_cast<int>(ExportList::DIFFU_ALPHA_X)]);
this->alpha_y = Rcpp::as<std::vector<double>>(
setup[static_cast<int>(ExportList::DIFFU_ALPHA_Y)]);
this->alpha_x =
Rcpp::List(setup[static_cast<int>(ExportList::DIFFU_ALPHA_X)]);
this->alpha_y =
Rcpp::List(setup[static_cast<int>(ExportList::DIFFU_ALPHA_Y)]);
this->database =
Rcpp::as<std::string>(setup[static_cast<int>(ExportList::CHEM_DATABASE)]);

13
src/Init/InitialList.hpp
View File

@ -1,5 +1,7 @@
#pragma once
#include <tug/Boundary.hpp>
#include <RInside.h>
#include <Rcpp.h>
#include <Rcpp/vector/instantiation.h>
@ -8,7 +10,9 @@
#include <cstdint>
#include <IPhreeqcPOET.hpp>
#include <map>
#include <string>
#include <utility>
#include <vector>
#include <DataStructures/Field.hpp>
@ -26,6 +30,8 @@ public:
void importList(const Rcpp::List &setup);
Rcpp::List exportList();
Field getInitialGrid() const { return Field(this->initial_grid); }
private:
RInside &R;
@ -93,6 +99,9 @@ private:
public:
struct DiffusionInit {
using BoundaryMap =
std::map<std::string, std::vector<tug::BoundaryElement<TugType>>>;
uint8_t dim;
std::uint32_t n_cols;
std::uint32_t n_rows;
@ -104,8 +113,8 @@ public:
std::vector<std::string> transport_names;
Field initial_grid;
Field boundaries;
BoundaryMap boundaries;
Field alpha_x;
Field alpha_y;
};

226
src/Transport/DiffusionModule.cpp
View File

@ -4,6 +4,9 @@
**
** Copyright (C) 2018-2022 Marco De Lucia, Max Luebke (GFZ Potsdam)
**
** Copyright (C) 2023-2024 Marco De Lucia (GFZ Potsdam), Max Luebke (University
** of Potsdam)
**
** POET is free software; you can redistribute it and/or modify it under the
** terms of the GNU General Public License as published by the Free Software
** Foundation; either version 2 of the License, or (at your option) any later
@ -19,174 +22,101 @@
*/
#include "DiffusionModule.hpp"
#include "../Base/Macros.hpp"
#include <tug/BoundaryCondition.hpp>
#include <tug/Diffusion.hpp>
#include "Base/Macros.hpp"
#include "Init/InitialList.hpp"
#include <Rcpp.h>
#include <algorithm>
#include <cstdint>
#include <Eigen/src/Core/Map.h>
#include <Eigen/src/Core/Matrix.h>
#include <Eigen/src/Core/util/Constants.h>
#include <chrono>
#include <array>
#include <cassert>
#include <mpi.h>
#include "tug/Boundary.hpp"
#include "tug/Grid.hpp"
#include "tug/Simulation.hpp"
#include <iostream>
#include <ostream>
#include <string>
#include <vector>
using namespace poet;
static constexpr double ZERO_MULTIPLICATOR = 10E-14;
static inline std::vector<TugType>
MatrixToVec(const Eigen::MatrixX<TugType> &mat) {
std::vector<TugType> vec(mat.rows() * mat.cols());
constexpr std::array<uint8_t, 4> borders = {
tug::bc::BC_SIDE_LEFT, tug::bc::BC_SIDE_RIGHT, tug::bc::BC_SIDE_TOP,
tug::bc::BC_SIDE_BOTTOM};
inline const char *convert_bc_to_config_file(uint8_t in) {
switch (in) {
case tug::bc::BC_SIDE_TOP:
return "N";
case tug::bc::BC_SIDE_RIGHT:
return "E";
case tug::bc::BC_SIDE_BOTTOM:
return "S";
case tug::bc::BC_SIDE_LEFT:
return "W";
}
return "";
}
DiffusionModule::DiffusionModule(const poet::DiffusionParams &diffu_args,
const poet::GridParams &grid_params)
: n_cells_per_prop(grid_params.total_n) {
this->diff_input.setGridCellN(grid_params.n_cells[0], grid_params.n_cells[1]);
this->diff_input.setDomainSize(grid_params.s_cells[0],
grid_params.s_cells[1]);
this->dim = grid_params.dim;
this->initialize(diffu_args, grid_params.total_n);
}
void DiffusionModule::initialize(const poet::DiffusionParams &args,
std::uint32_t n_grid_cells) {
// const poet::DiffusionParams args(this->R);
// name of props
//
this->prop_names = Rcpp::as<std::vector<std::string>>(args.initial_t.names());
this->prop_count = this->prop_names.size();
// initialize field and alphas
this->alpha.reserve(this->prop_count);
std::vector<std::vector<double>> initial_values;
for (uint32_t i = 0; i < this->prop_count; i++) {
// get alphas - we cannot assume alphas are provided in same order as
// initial input
this->alpha.push_back(args.alpha[this->prop_names[i]]);
const double val = args.initial_t[prop_names[i]];
std::vector<double> init_val(n_grid_cells, val);
initial_values.push_back(std::move(init_val));
if (this->dim == this->DIM_2D) {
tug::bc::BoundaryCondition bc(diff_input.grid.grid_cells[0],
diff_input.grid.grid_cells[1]);
this->bc_vec.push_back(bc);
} else {
tug::bc::BoundaryCondition bc(diff_input.grid.grid_cells[0]);
this->bc_vec.push_back(bc);
for (std::uint32_t i = 0; i < mat.cols(); i++) {
for (std::uint32_t j = 0; j < mat.rows(); j++) {
vec[j * mat.cols() + i] = mat(j, i);
}
}
this->t_field = Field(n_grid_cells, initial_values, prop_names);
// apply boundary conditions to each ghost node
uint8_t bc_size = (this->dim == this->DIM_1D ? 2 : 4);
for (uint8_t i = 0; i < bc_size; i++) {
const auto &side = borders[i];
std::vector<uint32_t> vecinj_i = Rcpp::as<std::vector<uint32_t>>(
args.vecinj_index[convert_bc_to_config_file(side)]);
for (int i = 0; i < this->prop_count; i++) {
std::vector<double> bc_vec = args.vecinj[this->prop_names[i]];
tug::bc::BoundaryCondition &curr_bc = *(this->bc_vec.begin() + i);
for (int j = 0; j < vecinj_i.size(); j++) {
if (vecinj_i[j] == 0) {
continue;
}
if (this->dim == this->DIM_2D) {
curr_bc(side, j) = {tug::bc::BC_TYPE_CONSTANT,
bc_vec[vecinj_i[j] - 1]};
}
if (this->dim == this->DIM_1D) {
curr_bc(side) = {tug::bc::BC_TYPE_CONSTANT, bc_vec[vecinj_i[j] - 1]};
}
}
}
}
// apply inner grid constant cells
// NOTE: opening a scope here for distinguish variable names
if (args.vecinj_inner.size() != 0) {
// apply inner constant cells for every concentration
for (int i = 0; i < this->prop_count; i++) {
std::vector<double> bc_vec = args.vecinj[this->prop_names[i]];
tug::bc::BoundaryCondition &curr_bc = *(this->bc_vec.begin() + i);
for (int j = 0; j < args.vecinj_inner.size(); j++) {
std::vector<double> inner_tuple =
Rcpp::as<std::vector<double>>(args.vecinj_inner[j]);
tug::bc::boundary_condition bc = {tug::bc::BC_TYPE_CONSTANT,
bc_vec[inner_tuple[0] - 1]};
this->index_constant_cells.push_back(inner_tuple[1]);
uint32_t x = inner_tuple[1];
uint32_t y = (this->dim == this->DIM_1D ? 0 : inner_tuple[2]);
curr_bc.setInnerBC(bc, x, y);
}
}
}
return vec;
}
void DiffusionModule::simulate(double dt) {
double sim_a_transport, sim_b_transport;
static inline Eigen::MatrixX<TugType>
VecToMatrix(const std::vector<TugType> &vec, std::uint32_t n_rows,
std::uint32_t n_cols) {
Eigen::MatrixX<TugType> mat(n_rows, n_cols);
sim_b_transport = MPI_Wtime();
MSG_NOENDL("DiffusionModule::simulate(): Starting diffusion ...");
std::cout << std::flush;
std::vector<std::vector<double>> field_2d = t_field.As2DVector();
this->diff_input.setTimestep(dt);
for (int i = 0; i < field_2d.size(); i++) {
std::vector<double> in_alpha(this->n_cells_per_prop, this->alpha[i]);
this->diff_input.setBoundaryCondition(this->bc_vec[i]);
if (this->dim == this->DIM_1D) {
tug::diffusion::BTCS_1D(this->diff_input, field_2d[i].data(),
in_alpha.data());
} else {
tug::diffusion::ADI_2D(this->diff_input, field_2d[i].data(),
in_alpha.data());
for (std::uint32_t i = 0; i < n_cols; i++) {
for (std::uint32_t j = 0; j < n_rows; j++) {
mat(j, i) = vec[j * n_cols + i];
}
}
t_field = field_2d;
sim_a_transport = MPI_Wtime();
transport_t += sim_a_transport - sim_b_transport;
std::cout << " done in " << sim_a_transport - sim_b_transport << "sec"
<< std::endl;
return mat;
}
void DiffusionModule::end() {
// R["simtime_transport"] = transport_t;
// R.parseEvalQ("profiling$simtime_transport <- simtime_transport");
// static constexpr double ZERO_MULTIPLICATOR = 10E-14;
void DiffusionModule::simulate(double requested_dt) {
MSG("Starting diffusion ...");
const auto start_diffusion_t = std::chrono::high_resolution_clock::now();
const auto &n_rows = this->param_list.n_rows;
const auto &n_cols = this->param_list.n_cols;
tug::Grid<TugType> grid(param_list.n_rows, param_list.n_cols);
tug::Boundary boundary(grid);
grid.setDomain(param_list.s_rows, param_list.s_cols);
tug::Simulation<TugType> sim(grid, boundary);
sim.setIterations(1);
for (const auto &sol_name : this->param_list.transport_names) {
auto &species_conc = this->transport_field[sol_name];
Eigen::MatrixX<TugType> conc = VecToMatrix(species_conc, n_rows, n_cols);
Eigen::MatrixX<TugType> alpha_x =
VecToMatrix(this->param_list.alpha_x[sol_name], n_rows, n_cols);
Eigen::MatrixX<TugType> alpha_y =
VecToMatrix(this->param_list.alpha_y[sol_name], n_rows, n_cols);
boundary.deserialize(this->param_list.boundaries[sol_name]);
grid.setAlpha(alpha_x, alpha_y);
grid.setConcentrations(conc);
sim.setTimestep(requested_dt);
sim.run();
species_conc = MatrixToVec(grid.getConcentrations());
}
const auto end_diffusion_t = std::chrono::high_resolution_clock::now();
const auto consumed_time_seconds =
std::chrono::duration_cast<std::chrono::duration<double>>(
end_diffusion_t - start_diffusion_t);
MSG("Diffusion done in " + std::to_string(consumed_time_seconds.count()) +
" [s]");
transport_t += consumed_time_seconds.count();
}
double DiffusionModule::getTransportTime() { return this->transport_t; }

49
src/Transport/DiffusionModule.hpp
View File

@ -21,20 +21,13 @@
#ifndef DIFFUSION_MODULE_H
#define DIFFUSION_MODULE_H
#include "../Base/Grid.hpp"
#include "../Base/SimParams.hpp"
#include "../DataStructures/DataStructures.hpp"
#include "DataStructures/Field.hpp"
#include "Init/InitialList.hpp"
#include <tug/BoundaryCondition.hpp>
#include <tug/Diffusion.hpp>
#include <array>
#include <cmath>
#include <cstdint>
#include <string>
#include <vector>
#include <sys/types.h>
namespace poet {
/**
* @brief Class describing transport simulation
*
@ -53,8 +46,8 @@ public:
*
* @param R RRuntime object
*/
DiffusionModule(const poet::DiffusionParams &diffu_args,
const poet::GridParams &grid_params);
DiffusionModule(const InitialList::DiffusionInit &init_list, Field field)
: param_list(init_list), transport_field(field){};
/**
* @brief Run simulation for one iteration
@ -64,14 +57,6 @@ public:
*/
void simulate(double dt);
/**
* @brief End simulation
*
* All measured timings are distributed to the R runtime
*
*/
void end();
/**
* @brief Get the transport time
*
@ -84,33 +69,17 @@ public:
*
* \return Reference to the diffusion field.
*/
Field &getField() { return this->t_field; }
Field &getField() { return this->transport_field; }
private:
/**
* @brief Instance of RRuntime
*
*/
// RRuntime &R;
enum { DIM_1D = 1, DIM_2D };
InitialList::DiffusionInit param_list;
void initialize(const poet::DiffusionParams &args,
std::uint32_t n_grid_cells);
uint8_t dim;
uint32_t prop_count;
tug::diffusion::TugInput diff_input;
std::vector<double> alpha;
std::vector<uint32_t> index_constant_cells;
std::vector<std::string> prop_names;
std::vector<tug::bc::BoundaryCondition> bc_vec;
Field t_field;
uint32_t n_cells_per_prop;
Field transport_field;
/**
* @brief time spent for transport

527
src/poet.cpp
View File

@ -4,6 +4,8 @@
**
** Copyright (C) 2018-2022 Marco De Lucia, Max Luebke (GFZ Potsdam)
**
** Copyright (C) 2023-2024 Max Luebke (University of Potsdam)
**
** POET is free software; you can redistribute it and/or modify it under the
** terms of the GNU General Public License as published by the Free Software
** Foundation; either version 2 of the License, or (at your option) any later
@ -18,125 +20,308 @@
** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
#include "Base/Grid.hpp"
#include "Base/Macros.hpp"
#include "Base/RInsidePOET.hpp"
#include "Base/SimParams.hpp"
#include "Chemistry/ChemistryModule.hpp"
// #include "Chemistry/ChemistryModule.hpp"
#include "DataStructures/Field.hpp"
#include "Init/InitialList.hpp"
#include "Transport/DiffusionModule.hpp"
#include <RInside.h>
#include <R_ext/Boolean.h>
#include <Rcpp/DataFrame.h>
#include <Rcpp/vector/instantiation.h>
#include <poet.hpp>
#include <Rcpp.h>
#include <cstdint>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <string>
#include <vector>
#include <mpi.h>
#include "Base/argh.hpp"
using namespace std;
using namespace poet;
using namespace Rcpp;
poet::ChemistryModule::SingleCMap DFToHashMap(const Rcpp::DataFrame &df) {
std::unordered_map<std::string, double> out_map;
vector<string> col_names = Rcpp::as<vector<string>>(df.names());
static int MY_RANK = 0;
for (const auto &name : col_names) {
double val = df[name.c_str()];
out_map.insert({name, val});
}
// poet::ChemistryModule::SingleCMap DFToHashMap(const Rcpp::DataFrame &df) {
// std::unordered_map<std::string, double> out_map;
// vector<string> col_names = Rcpp::as<vector<string>>(df.names());
return out_map;
}
// for (const auto &name : col_names) {
// double val = df[name.c_str()];
// out_map.insert({name, val});
// }
// return out_map;
// }
// 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) {
R["TMP"] = Rcpp::wrap(trans.AsVector());
R["TMP_PROPS"] = Rcpp::wrap(trans.GetProps());
R.parseEval(std::string(
"mysetup$state_T <- setNames(data.frame(matrix(TMP, nrow=" +
std::to_string(trans.GetRequestedVecSize()) + ")), TMP_PROPS)"));
void writeFieldsToR(RInside &R, const Field &trans) {
// , const Field &chem) {
R["TMP"] = Rcpp::wrap(chem.AsVector());
R["TMP_PROPS"] = Rcpp::wrap(chem.GetProps());
R.parseEval(std::string(
"mysetup$state_C <- setNames(data.frame(matrix(TMP, nrow=" +
std::to_string(chem.GetRequestedVecSize()) + ")), TMP_PROPS)"));
Rcpp::DataFrame t_field(trans.asSEXP());
R["TMP"] = t_field;
R.parseEval("mysetup$state_T <- TMP");
R.parseEval("mysetup$state_C <- TMP");
// R["TMP"] = Rcpp::wrap(trans.AsVector());
// R["TMP_PROPS"] = Rcpp::wrap(trans.GetProps());
// R.parseEval(std::string(
// "mysetup$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(
// "mysetup$state_C <- setNames(data.frame(matrix(TMP, nrow=" +
// std::to_string(chem.GetRequestedVecSize()) + ")), TMP_PROPS)"));
}
void set_chem_parameters(poet::ChemistryModule &chem, uint32_t wp_size,
const std::string &database_path) {
chem.SetErrorHandlerMode(1);
chem.SetComponentH2O(false);
chem.SetRebalanceFraction(0.5);
chem.SetRebalanceByCell(true);
chem.UseSolutionDensityVolume(false);
chem.SetPartitionUZSolids(false);
enum ParseRet { PARSER_OK, PARSER_ERROR, PARSER_HELP };
// Set concentration units
// 1, mg/L; 2, mol/L; 3, kg/kgs
chem.SetUnitsSolution(2);
// 0, mol/L cell; 1, mol/L water; 2 mol/L rock
chem.SetUnitsPPassemblage(1);
// 0, mol/L cell; 1, mol/L water; 2 mol/L rock
chem.SetUnitsExchange(1);
// 0, mol/L cell; 1, mol/L water; 2 mol/L rock
chem.SetUnitsSurface(1);
// 0, mol/L cell; 1, mol/L water; 2 mol/L rock
chem.SetUnitsGasPhase(1);
// 0, mol/L cell; 1, mol/L water; 2 mol/L rock
chem.SetUnitsSSassemblage(1);
// 0, mol/L cell; 1, mol/L water; 2 mol/L rock
chem.SetUnitsKinetics(1);
ParseRet parseInitValues(char **argv, RInsidePOET &R,
RuntimeParameters &params) {
// initialize argh object
argh::parser cmdl(argv);
// Set representative volume
std::vector<double> rv;
rv.resize(wp_size, 1.0);
chem.SetRepresentativeVolume(rv);
// if user asked for help
if (cmdl[{"help", "h"}]) {
if (MY_RANK == 0) {
MSG("Todo");
MSG("See README.md for further information.");
}
// Set initial porosity
std::vector<double> por;
por.resize(wp_size, 1);
chem.SetPorosity(por);
return ParseRet::PARSER_HELP;
}
// if positional arguments are missing
else if (!cmdl(2)) {
if (MY_RANK == 0) {
ERRMSG("POET needs 2 positional arguments: ");
ERRMSG("1) the R script defining your simulation and");
ERRMSG("2) the directory prefix where to save results and profiling");
}
return ParseRet::PARSER_ERROR;
}
// Set initial saturation
std::vector<double> sat;
sat.resize(wp_size, 1.0);
chem.SetSaturation(sat);
// parse flags and check for unknown
for (const auto &option : cmdl.flags()) {
if (!(flaglist.find(option) != flaglist.end())) {
if (MY_RANK == 0) {
ERRMSG("Unrecognized option: " + option);
ERRMSG("Make sure to use available options. Exiting!");
}
return ParseRet::PARSER_ERROR;
}
}
// Load database
chem.LoadDatabase(database_path);
// parse parameters and check for unknown
for (const auto &option : cmdl.params()) {
if (!(paramlist.find(option.first) != paramlist.end())) {
if (MY_RANK == 0) {
ERRMSG("Unrecognized option: " + option.first);
ERRMSG("Make sure to use available options. Exiting!");
}
return ParseRet::PARSER_ERROR;
}
}
// simparams.print_progressbar = cmdl[{"P", "progress"}];
// // simparams.print_progressbar = cmdl[{"P", "progress"}];
// /* Parse DHT arguments */
// chem_params.use_dht = cmdl["dht"];
// chem_params.use_interp = cmdl["interp"];
// // cout << "CPP: DHT is " << ( dht_enabled ? "ON" : "OFF" ) << '\n';
// cmdl("dht-size", DHT_SIZE_PER_PROCESS_MB) >> chem_params.dht_size;
// // cout << "CPP: DHT size per process (Byte) = " << dht_size_per_process <<
// // endl;
// cmdl("dht-snaps", 0) >> chem_params.dht_snaps;
// cmdl("dht-file") >> chem_params.dht_file;
// /*Parse work package size*/
// cmdl("work-package-size", WORK_PACKAGE_SIZE_DEFAULT) >> simparams.wp_size;
// cmdl("interp-size", 100) >> chem_params.pht_size;
// cmdl("interp-min", 5) >> chem_params.interp_min_entries;
// cmdl("interp-bucket-entries", 20) >> chem_params.pht_max_entries;
// /*Parse output options*/
// simparams.store_result = !cmdl["ignore-result"];
// /*Parse work package size*/
// cmdl("work-package-size", WORK_PACKAGE_SIZE_DEFAULT) >> simparams.wp_size;
// chem_params.use_interp = cmdl["interp"];
// cmdl("interp-size", 100) >> chem_params.pht_size;
// cmdl("interp-min", 5) >> chem_params.interp_min_entries;
// cmdl("interp-bucket-entries", 20) >> chem_params.pht_max_entries;
// /*Parse output options*/
// simparams.store_result = !cmdl["ignore-result"];
// if (simparams.world_rank == 0) {
// MSG("Complete results storage is " + BOOL_PRINT(simparams.store_result));
// MSG("Work Package Size: " + std::to_string(simparams.wp_size));
// MSG("DHT is " + BOOL_PRINT(chem_params.use_dht));
// if (chem_params.use_dht) {
// MSG("DHT strategy is " + std::to_string(simparams.dht_strategy));
// // MDL: these should be outdated (?)
// // MSG("DHT key default digits (ignored if 'signif_vector' is "
// // "defined) = "
// // << simparams.dht_significant_digits);
// // MSG("DHT logarithm before rounding: "
// // << (simparams.dht_log ? "ON" : "OFF"));
// MSG("DHT size per process (Megabyte) = " +
// std::to_string(chem_params.dht_size));
// MSG("DHT save snapshots is " + BOOL_PRINT(chem_params.dht_snaps));
// MSG("DHT load file is " + chem_params.dht_file);
// }
// if (chem_params.use_interp) {
// MSG("PHT interpolation enabled: " +
// BOOL_PRINT(chem_params.use_interp)); MSG("PHT interp-size = " +
// std::to_string(chem_params.pht_size)); MSG("PHT interp-min = " +
// std::to_string(chem_params.interp_min_entries));
// MSG("PHT interp-bucket-entries = " +
// std::to_string(chem_params.pht_max_entries));
// }
// }
std::string init_file;
std::string runtime_file;
std::string out_dir;
cmdl(1) >> init_file;
cmdl(2) >> runtime_file;
cmdl(3) >> out_dir;
// chem_params.dht_outdir = out_dir;
/* distribute information to R runtime */
// if local_rank == 0 then master else worker
R["local_rank"] = MY_RANK;
// assign a char* (string) to 'filesim'
R["filesim"] = wrap(runtime_file);
// assign a char* (string) to 'fileout'
R["fileout"] = wrap(out_dir);
// pass the boolean "store_result" to the R process
// R["store_result"] = simparams.store_result;
// // worker count
// R["n_procs"] = simparams.world_size - 1;
// // work package size
// R["work_package_size"] = simparams.wp_size;
// // dht enabled?
// R["dht_enabled"] = chem_params.use_dht;
// // log before rounding?
// R["dht_log"] = simparams.dht_log;
try {
Rcpp::Function source("source");
Rcpp::Function readRDS("readRDS");
Rcpp::List init_params_ = readRDS(init_file);
params.init_params = init_params_;
Rcpp::List runtime_params =
source(runtime_file, Rcpp::Named("local", true));
runtime_params = runtime_params["value"];
R[r_runtime_parameters] = runtime_params;
params.timesteps =
Rcpp::as<std::vector<double>>(runtime_params["timesteps"]);
} catch (const std::exception &e) {
ERRMSG("Error while parsing R scripts: " + std::string(e.what()));
return ParseRet::PARSER_ERROR;
}
// eval the init string, ignoring any returns
// R.parseEvalQ("source(filesim)");
// R.parseEvalQ("mysetup <- setup");
// this->chem_params.initFromR(R);
return ParseRet::PARSER_OK;
}
inline double RunMasterLoop(SimParams &params, RInside &R,
const GridParams &g_params, uint32_t nxyz_master) {
// void set_chem_parameters(poet::ChemistryModule &chem, uint32_t wp_size,
// const std::string &database_path) {
// chem.SetErrorHandlerMode(1);
// chem.SetComponentH2O(false);
// chem.SetRebalanceFraction(0.5);
// chem.SetRebalanceByCell(true);
// chem.UseSolutionDensityVolume(false);
// chem.SetPartitionUZSolids(false);
// // Set concentration units
// // 1, mg/L; 2, mol/L; 3, kg/kgs
// chem.SetUnitsSolution(2);
// // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
// chem.SetUnitsPPassemblage(1);
// // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
// chem.SetUnitsExchange(1);
// // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
// chem.SetUnitsSurface(1);
// // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
// chem.SetUnitsGasPhase(1);
// // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
// chem.SetUnitsSSassemblage(1);
// // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
// chem.SetUnitsKinetics(1);
// // Set representative volume
// std::vector<double> rv;
// rv.resize(wp_size, 1.0);
// chem.SetRepresentativeVolume(rv);
// // Set initial porosity
// std::vector<double> por;
// por.resize(wp_size, 1);
// chem.SetPorosity(por);
// // Set initial saturation
// std::vector<double> sat;
// sat.resize(wp_size, 1.0);
// chem.SetSaturation(sat);
// // Load database
// chem.LoadDatabase(database_path);
// }
static double RunMasterLoop(RInside &R, const RuntimeParameters &params,
DiffusionModule &diffusion) {
DiffusionParams d_params{R};
DiffusionModule diffusion(d_params, g_params);
/* Iteration Count is dynamic, retrieving value from R (is only needed by
* master for the following loop) */
uint32_t maxiter = R.parseEval("mysetup$iterations");
uint32_t maxiter = params.timesteps.size();
double sim_time = .0;
ChemistryModule chem(nxyz_master, params.getNumParams().wp_size, maxiter,
params.getChemParams(), MPI_COMM_WORLD);
// ChemistryModule chem(nxyz_master, params.getNumParams().wp_size, maxiter,
// params.getChemParams(), MPI_COMM_WORLD);
set_chem_parameters(chem, nxyz_master, params.getChemParams().database_path);
chem.RunInitFile(params.getChemParams().input_script);
// set_chem_parameters(chem, nxyz_master,
// params.getChemParams().database_path);
// chem.RunInitFile(params.getChemParams().input_script);
poet::ChemistryModule::SingleCMap init_df = DFToHashMap(d_params.initial_t);
chem.initializeField(diffusion.getField());
// poet::ChemistryModule::SingleCMap init_df =
// DFToHashMap(d_params.initial_t);
// chem.initializeField(diffusion.getField());
if (params.getNumParams().print_progressbar) {
chem.setProgressBarPrintout(true);
}
// if (params.getNumParams().print_progressbar) {
// chem.setProgressBarPrintout(true);
// }
/* SIMULATION LOOP */
@ -147,31 +332,30 @@ inline double RunMasterLoop(SimParams &params, RInside &R,
// cout << "CPP: Evaluating next time step" << endl;
// R.parseEvalQ("mysetup <- master_iteration_setup(mysetup)");
double dt = Rcpp::as<double>(
R.parseEval("mysetup$timesteps[" + std::to_string(iter) + "]"));
double dt = params.timesteps[iter - 1];
// cout << "CPP: Next time step is " << dt << "[s]" << endl;
MSG("Next time step is " + std::to_string(dt) + " [s]");
/* displaying iteration number, with C++ and R iterator */
MSG("Going through iteration " + std::to_string(iter));
MSG("R's $iter: " +
std::to_string((uint32_t)(R.parseEval("mysetup$iter"))) +
". Iteration");
// MSG("R's $iter: " +
// std::to_string((uint32_t)(R.parseEval("mysetup$iter"))) +
// ". Iteration");
/* run transport */
// TODO: transport to diffusion
diffusion.simulate(dt);
chem.getField().update(diffusion.getField());
// chem.getField().update(diffusion.getField());
MSG("Chemistry step");
chem.SetTimeStep(dt);
chem.RunCells();
// chem.SetTimeStep(dt);
// chem.RunCells();
writeFieldsToR(R, diffusion.getField(), chem.GetField());
diffusion.getField().update(chem.GetField());
writeFieldsToR(R, diffusion.getField());
// diffusion.getField().update(chem.GetField());
R["req_dt"] = dt;
R["simtime"] = (sim_time += dt);
@ -183,7 +367,7 @@ inline double RunMasterLoop(SimParams &params, RInside &R,
// state_C after every iteration if the cmdline option
// --ignore-results is not given (and thus the R variable
// store_result is TRUE)
R.parseEvalQ("mysetup <- master_iteration_end(setup=mysetup)");
R.parseEval("mysetup <- master_iteration_end(setup=mysetup)");
MSG("End of *coupling* iteration " + std::to_string(iter) + "/" +
std::to_string(maxiter));
@ -194,60 +378,60 @@ inline double RunMasterLoop(SimParams &params, RInside &R,
dSimTime += end_t - start_t;
} // END SIMULATION LOOP
R.parseEvalQ("profiling <- list()");
// R.parseEvalQ("profiling <- list()");
R["simtime_chemistry"] = chem.GetChemistryTime();
R.parseEvalQ("profiling$simtime_chemistry <- simtime_chemistry");
// R["simtime_chemistry"] = chem.GetChemistryTime();
// R.parseEvalQ("profiling$simtime_chemistry <- simtime_chemistry");
R["chemistry_loop"] = chem.GetMasterLoopTime();
R.parseEvalQ("profiling$chemistry_loop <- chemistry_loop");
// R["chemistry_loop"] = chem.GetMasterLoopTime();
// R.parseEvalQ("profiling$chemistry_loop <- chemistry_loop");
R["chemistry_sequential"] = chem.GetMasterSequentialTime();
R.parseEvalQ("profiling$simtime_sequential <- chemistry_sequential");
// R["chemistry_sequential"] = chem.GetMasterSequentialTime();
// R.parseEvalQ("profiling$simtime_sequential <- chemistry_sequential");
R["idle_master"] = chem.GetMasterIdleTime();
R.parseEvalQ("profiling$idle_master <- idle_master");
// R["idle_master"] = chem.GetMasterIdleTime();
// R.parseEvalQ("profiling$idle_master <- idle_master");
R["idle_worker"] = Rcpp::wrap(chem.GetWorkerIdleTimings());
R.parseEvalQ("profiling$idle_worker <- idle_worker");
// R["idle_worker"] = Rcpp::wrap(chem.GetWorkerIdleTimings());
// R.parseEvalQ("profiling$idle_worker <- idle_worker");
R["phreeqc_time"] = Rcpp::wrap(chem.GetWorkerPhreeqcTimings());
R.parseEvalQ("profiling$phreeqc <- phreeqc_time");
// R["phreeqc_time"] = Rcpp::wrap(chem.GetWorkerPhreeqcTimings());
// R.parseEvalQ("profiling$phreeqc <- phreeqc_time");
R["simtime_transport"] = diffusion.getTransportTime();
R.parseEvalQ("profiling$simtime_transport <- simtime_transport");
// R["simtime_transport"] = diffusion.getTransportTime();
// R.parseEvalQ("profiling$simtime_transport <- simtime_transport");
// R["phreeqc_count"] = phreeqc_counts;
// R.parseEvalQ("profiling$phreeqc_count <- phreeqc_count");
if (params.getChemParams().use_dht) {
R["dht_hits"] = Rcpp::wrap(chem.GetWorkerDHTHits());
R.parseEvalQ("profiling$dht_hits <- dht_hits");
R["dht_evictions"] = Rcpp::wrap(chem.GetWorkerDHTEvictions());
R.parseEvalQ("profiling$dht_evictions <- dht_evictions");
R["dht_get_time"] = Rcpp::wrap(chem.GetWorkerDHTGetTimings());
R.parseEvalQ("profiling$dht_get_time <- dht_get_time");
R["dht_fill_time"] = Rcpp::wrap(chem.GetWorkerDHTFillTimings());
R.parseEvalQ("profiling$dht_fill_time <- dht_fill_time");
}
if (params.getChemParams().use_interp) {
R["interp_w"] = Rcpp::wrap(chem.GetWorkerInterpolationWriteTimings());
R.parseEvalQ("profiling$interp_write <- interp_w");
R["interp_r"] = Rcpp::wrap(chem.GetWorkerInterpolationReadTimings());
R.parseEvalQ("profiling$interp_read <- interp_r");
R["interp_g"] = Rcpp::wrap(chem.GetWorkerInterpolationGatherTimings());
R.parseEvalQ("profiling$interp_gather <- interp_g");
R["interp_fc"] =
Rcpp::wrap(chem.GetWorkerInterpolationFunctionCallTimings());
R.parseEvalQ("profiling$interp_function_calls <- interp_fc");
R["interp_calls"] = Rcpp::wrap(chem.GetWorkerInterpolationCalls());
R.parseEvalQ("profiling$interp_calls <- interp_calls");
R["interp_cached"] = Rcpp::wrap(chem.GetWorkerPHTCacheHits());
R.parseEvalQ("profiling$interp_cached <- interp_cached");
}
// if (params.getChemParams().use_dht) {
// R["dht_hits"] = Rcpp::wrap(chem.GetWorkerDHTHits());
// R.parseEvalQ("profiling$dht_hits <- dht_hits");
// R["dht_evictions"] = Rcpp::wrap(chem.GetWorkerDHTEvictions());
// R.parseEvalQ("profiling$dht_evictions <- dht_evictions");
// R["dht_get_time"] = Rcpp::wrap(chem.GetWorkerDHTGetTimings());
// R.parseEvalQ("profiling$dht_get_time <- dht_get_time");
// R["dht_fill_time"] = Rcpp::wrap(chem.GetWorkerDHTFillTimings());
// R.parseEvalQ("profiling$dht_fill_time <- dht_fill_time");
// }
// if (params.getChemParams().use_interp) {
// R["interp_w"] = Rcpp::wrap(chem.GetWorkerInterpolationWriteTimings());
// R.parseEvalQ("profiling$interp_write <- interp_w");
// R["interp_r"] = Rcpp::wrap(chem.GetWorkerInterpolationReadTimings());
// R.parseEvalQ("profiling$interp_read <- interp_r");
// R["interp_g"] =
// Rcpp::wrap(chem.GetWorkerInterpolationGatherTimings());
// R.parseEvalQ("profiling$interp_gather <- interp_g");
// R["interp_fc"] =
// Rcpp::wrap(chem.GetWorkerInterpolationFunctionCallTimings());
// R.parseEvalQ("profiling$interp_function_calls <- interp_fc");
// R["interp_calls"] = Rcpp::wrap(chem.GetWorkerInterpolationCalls());
// R.parseEvalQ("profiling$interp_calls <- interp_calls");
// R["interp_cached"] = Rcpp::wrap(chem.GetWorkerPHTCacheHits());
// R.parseEvalQ("profiling$interp_cached <- interp_cached");
// }
chem.MasterLoopBreak();
diffusion.end();
// chem.MasterLoopBreak();
return dSimTime;
}
@ -271,24 +455,24 @@ int main(int argc, char *argv[]) {
if (world_rank > 0) {
{
SimParams params(world_rank, world_size);
int pret = params.parseFromCmdl(argv, R);
// SimParams params(world_rank, world_size);
// int pret = params.parseFromCmdl(argv, R);
if (pret == poet::PARSER_ERROR) {
MPI_Finalize();
return EXIT_FAILURE;
} else if (pret == poet::PARSER_HELP) {
MPI_Finalize();
return EXIT_SUCCESS;
}
// if (pret == poet::PARSER_ERROR) {
// MPI_Finalize();
// return EXIT_FAILURE;
// } else if (pret == poet::PARSER_HELP) {
// MPI_Finalize();
// return EXIT_SUCCESS;
// }
// ChemistryModule worker(nxyz, nxyz, MPI_COMM_WORLD);
ChemistryModule worker = poet::ChemistryModule::createWorker(
MPI_COMM_WORLD, params.getChemParams());
set_chem_parameters(worker, worker.GetWPSize(),
params.getChemParams().database_path);
// // ChemistryModule worker(nxyz, nxyz, MPI_COMM_WORLD);
// ChemistryModule worker = poet::ChemistryModule::createWorker(
// MPI_COMM_WORLD, params.getChemParams());
// set_chem_parameters(worker, worker.GetWPSize(),
// params.getChemParams().database_path);
worker.WorkerLoop();
// worker.WorkerLoop();
}
MPI_Barrier(MPI_COMM_WORLD);
@ -304,33 +488,31 @@ int main(int argc, char *argv[]) {
// TODO: kann raus
R.parseEvalQ(kin_r_library);
SimParams params(world_rank, world_size);
int pret = params.parseFromCmdl(argv, R);
RuntimeParameters run_params;
if (pret == poet::PARSER_ERROR) {
switch (parseInitValues(argv, R, run_params)) {
case ParseRet::PARSER_ERROR:
case ParseRet::PARSER_HELP:
MPI_Finalize();
return EXIT_FAILURE;
} else if (pret == poet::PARSER_HELP) {
MPI_Finalize();
return EXIT_SUCCESS;
return 0;
case ParseRet::PARSER_OK:
break;
}
MSG("RInside initialized on process " + std::to_string(world_rank));
R.parseEvalQ("mysetup <- setup");
// if (world_rank == 0) { // get timestep vector from
// grid_init function ... //
std::string master_init_code = "mysetup <- master_init(setup=setup)";
// R.parseEvalQ("mysetup <- setup");
// // if (world_rank == 0) { // get timestep vector from
// // grid_init function ... //
std::string master_init_code = "mysetup <- master_init(setup=mysetup)";
R.parseEval(master_init_code);
GridParams g_params(R);
params.initVectorParams(R);
// run_params.initVectorParams(R);
// MDL: store all parameters
if (world_rank == 0) {
MSG("Calling R Function to store calling parameters");
R.parseEvalQ("StoreSetup(setup=mysetup)");
// R.parseEvalQ("StoreSetup(setup=mysetup)");
}
if (world_rank == 0) {
@ -339,23 +521,28 @@ int main(int argc, char *argv[]) {
// MPI_Barrier(MPI_COMM_WORLD);
uint32_t nxyz_master = (world_size == 1 ? g_params.total_n : 1);
InitialList init_list(R);
init_list.importList(run_params.init_params);
dSimTime = RunMasterLoop(params, R, g_params, nxyz_master);
DiffusionModule diffusion(init_list.getDiffusionInit(),
init_list.getInitialGrid());
dSimTime = RunMasterLoop(R, run_params, diffusion);
MSG("finished simulation loop");
MSG("start timing profiling");
R["simtime"] = dSimTime;
R.parseEvalQ("profiling$simtime <- simtime");
// R["simtime"] = dSimTime;
// R.parseEvalQ("profiling$simtime <- simtime");
string r_vis_code;
r_vis_code = "saveRDS(profiling, file=paste0(fileout,'/timings.rds'));";
R.parseEval(r_vis_code);
// string r_vis_code;
// r_vis_code = "saveRDS(profiling, file=paste0(fileout,'/timings.rds'));";
// R.parseEval(r_vis_code);
MSG("Done! Results are stored as R objects into <" + params.getOutDir() +
"/timings.rds>");
// MSG("Done! Results are stored as R objects into <" + run_params.getOutDir()
// +
// "/timings.rds>");
MPI_Barrier(MPI_COMM_WORLD);

77
src/poet.hpp.in
View File

@ -1,12 +1,79 @@
#ifndef POET_H
#define POET_H
/*
** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
** Potsdam)
**
** Copyright (C) 2018-2023 Marco De Lucia, Max Luebke (GFZ Potsdam)
**
** Copyright (C) 2023-2024 Max Luebke (University of Potsdam)
**
** POET is free software; you can redistribute it and/or modify it under the
** terms of the GNU General Public License as published by the Free Software
** Foundation; either version 2 of the License, or (at your option) any later
** version.
**
** POET is distributed in the hope that it will be useful, but WITHOUT ANY
** WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
** A PARTICULAR PURPOSE. See the GNU General Public License for more details.
**
** You should have received a copy of the GNU General Public License along with
** this program; if not, write to the Free Software Foundation, Inc., 51
** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
#pragma once
#include <set>
#include <string>
#include <vector>
#include <Rcpp.h>
static const char *poet_version = "@POET_VERSION@";
// using the Raw string literal to avoid escaping the quotes
static inline std::string kin_r_library = R"(@R_KIN_LIB@)";
static const inline std::string kin_r_library = R"(@R_KIN_LIB@)";
static inline std::string init_r_library = R"(@R_INIT_LIB@)";
static const inline std::string init_r_library = R"(@R_INIT_LIB@)";
#endif // POET_H
static const inline std::string r_runtime_parameters = "mysetup";
const std::set<std::string> flaglist{"ignore-result", "dht", "P", "progress",
"interp"};
const std::set<std::string> paramlist{
"work-package-size", "dht-strategy", "dht-size", "dht-snaps",
"dht-file", "interp-size", "interp-min", "interp-bucket-entries"};
struct RuntimeParameters {
std::vector<double> timesteps;
bool print_progressbar;
Rcpp::List init_params;
struct ChemistryParams {
// std::string database_path;
// std::string input_script;
// bool use_dht;
// std::uint64_t dht_size;
// int dht_snaps;
// std::string dht_file;
// std::string dht_outdir;
// NamedVector<std::uint32_t> dht_signifs;
// bool use_interp;
// std::uint64_t pht_size;
// std::uint32_t pht_max_entries;
// std::uint32_t interp_min_entries;
// NamedVector<std::uint32_t> pht_signifs;
// struct Chem_Hook_Functions {
// RHookFunction<bool> dht_fill;
// RHookFunction<std::vector<double>> dht_fuzz;
// RHookFunction<std::vector<std::size_t>> interp_pre;
// RHookFunction<bool> interp_post;
// } hooks;
// void initFromR(RInsidePOET &R);
};
};