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

Merge branch 'ml/surfex' into 'main'

Browse Source 
Enable Surfex/Exchange reactants using IPhreeqc/POET API

See merge request naaice/poet!27
This commit is contained in:
Max Lübke 2024-05-06 12:11:50 +02:00
commit cc45bbc191
87 changed files with 3142 additions and 3664 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
.devcontainer/Dockerfile
View File

@ -2,11 +2,19 @@ FROM mcr.microsoft.com/vscode/devcontainers/base:debian
RUN sudo apt-get update && export DEBIAN_FRONTEND=noninteractive \
&& sudo apt-get install -y \
cmake-curses-gui \
clangd \
cmake \
git \
libeigen3-dev \
libopenmpi-dev \
r-cran-rcpp \
r-cran-rinside
r-base-dev
RUN git clone https://github.com/doctest/doctest.git /doctest \
&& cd /doctest \
&& mkdir build \
&& cd build \
&& cmake .. \
&& make install \
&& cd / \
&& rm -rf /doctest
RUN /usr/bin/R -q -e "install.packages(c('Rcpp', 'RInside'))"

3
.gitignore vendored
View File

@ -142,4 +142,5 @@ vignettes/*.pdf
build/
/.cache/
.vscode
.vscode
.codechecker

9
.gitmodules vendored
View File

@ -2,9 +2,6 @@
path = ext/tug
url = ../tug.git
[submodule "ext/phreeqcrm"]
path = ext/phreeqcrm
url = ../phreeqcrm-gfz.git
[submodule "ext/doctest"]
path = ext/doctest
url = https://github.com/doctest/doctest.git
[submodule "ext/iphreeqc"]
path = ext/iphreeqc
url = ../iphreeqc.git

12
CMakeLists.txt
View File

@ -1,12 +1,11 @@
# Version 3.9+ offers new MPI package variables
cmake_minimum_required(VERSION 3.9)
cmake_minimum_required(VERSION 3.14)
project(POET
LANGUAGES CXX C
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)
@ -16,8 +15,6 @@ list(APPEND CMAKE_MODULE_PATH "${POET_SOURCE_DIR}/CMake")
get_poet_version()
# set(GCC_CXX_FLAGS "-D STRICT_R_HEADERS") add_definitions(${GCC_CXX_FLAGS})
find_package(MPI REQUIRED)
find_package(RRuntime REQUIRED)
@ -29,15 +26,14 @@ add_subdirectory(bench)
set(TUG_ENABLE_TESTING OFF CACHE BOOL "" FORCE)
add_subdirectory(ext/tug EXCLUDE_FROM_ALL)
add_subdirectory(ext/phreeqcrm EXCLUDE_FROM_ALL)
add_subdirectory(ext/iphreeqc EXCLUDE_FROM_ALL)
option(POET_ENABLE_TESTING "Build test suite for POET" OFF)
if (POET_ENABLE_TESTING)
add_subdirectory(ext/doctest EXCLUDE_FROM_ALL)
add_subdirectory(test)
endif()
option(BUILD_DOC "Build documentation with doxygen" OFF)
add_subdirectory(docs)
add_subdirectory(docs)

108
README.md
View File

@ -20,10 +20,10 @@ pages](https://naaice.git-pages.gfz-potsdam.de/poet).
The following external header library is shipped with POET:
- **argh** - https://github.com/adishavit/argh (BSD license)
- **PhreeqcRM** with patches from GFZ -
https://www.usgs.gov/software/phreeqc-version-3 -
https://git.gfz-potsdam.de/mluebke/phreeqcrm-gfz
- **tug** - https://git.gfz-potsdam.de/sec34/tug
- **IPhreeqc** with patches from GFZ -
https://github.com/usgs-coupled/iphreeqc -
https://git.gfz-potsdam.de/naaice/iphreeqc
- **tug** - https://git.gfz-potsdam.de/naaice/tug
## Installation
@ -35,6 +35,7 @@ To compile POET you need several software to be installed:
- MPI-Implementation (tested with OpenMPI and MVAPICH)
- R language and environment
- CMake 3.9+
- Eigen3 3.4+ (required by `tug`)
- *optional*: `doxygen` with `dot` bindings for documentiation
The following R libraries must then be installed, which will get the
@ -107,58 +108,50 @@ The correspondending directory tree would look like this:
```sh
poet
├── bin
│ └── poet
├── R_lib
│ └── kin_r_library.R
│   ├── poet
│   └── poet_init
└── share
└── poet
└── bench
├── barite
│ ├── barite_interp_eval.R
│ ├── barite.pqi
│ ├── barite.R
│ └── db_barite.dat
├── dolo
│ ├── dolo_diffu_inner_large.R
│ ├── dolo_diffu_inner.R
│ ├── dolo_inner.pqi
│ ├── dolo_interp_long.R
│ └── phreeqc_kin.dat
└── surfex
├── ExBase.pqi
├── ex.R
├── SMILE_2021_11_01_TH.dat
├── SurfExBase.pqi
└── surfex.R
├── barite
│   ├── barite_200.rds
│   ├── barite_200_rt.R
│   ├── barite_het.rds
│   └── barite_het_rt.R
├── dolo
│   ├── dolo_inner_large.rds
│   ├── dolo_inner_large_rt.R
│   ├── dolo_interp.rds
│   └── dolo_interp_rt.R
└── surfex
├── PoetEGU_surfex_500.rds
└── PoetEGU_surfex_500_rt.R
```
The R libraries will be loaded at runtime and the paths are hardcoded
absolute paths inside `poet.cpp`. So, if you consider to move
`bin/poet` either change paths of the R source files and recompile
POET or also move `R_lib/*` relative to the binary.
With the installation of POET, two executables are provided:
- `poet` - the main executable to run simulations
- `poet_init` - a preprocessor to generate input files for POET from R scripts
The benchmarks consist of input scripts, which are provided as .R files.
Additionally, Phreeqc scripts and their corresponding databases are required,
stored as .pqi and .dat files, respectively.
Preprocessed benchmarks can be found in the `share/poet` directory with an
according *runtime* setup. More on those files and how to create them later.
## Running
Run POET by `mpirun ./poet <OPTIONS> <SIMFILE> <OUTPUT_DIRECTORY>`
Run POET by `mpirun ./poet [OPTIONS] <RUNFILE> <SIMFILE> <OUTPUT_DIRECTORY>`
where:
- **OPTIONS** - runtime parameters (explained below)
- **SIMFILE** - simulation described as R script (e.g.
`<POET_INSTALL_DIR>/share/poet/bench/dolo/dolo_interp_long.R`)
- **OPTIONS** - POET options (explained below)
- **RUNFILE** - Runtime parameters described as R script
- **SIMFILE** - Simulation input prepared by `poet_init`
- **OUTPUT_DIRECTORY** - path, where all output of POET should be stored
### Runtime options
### POET options
The following parameters can be set:
| Option | Value | Description |
|-----------------------------|--------------|--------------------------------------------------------------------------------------------------------------------------|
| **--work-package-size=** | _1..n_ | size of work packages (defaults to _5_) |
| **--ignore-result** | | disables store of simulation resuls |
| **-P, --progress** | | show progress bar |
| **--dht** | | enabling DHT usage (defaults to _OFF_) |
| **--dht-strategy=** | _0-1_ | change DHT strategy. **NOT IMPLEMENTED YET** (Defaults to _0_) |
| **--dht-size=** | _1-n_ | size of DHT per process involved in megabyte (defaults to _1000 MByte_) |
@ -180,14 +173,16 @@ Following values can be set:
### Example: Running from scratch
We will continue the above example and start a simulation with
`dolo_diffu_inner.R`. As transport a simple fixed-coefficient diffusion is used.
It's a 2D, 100x100 grid, simulating 10 time steps. To start the simulation with
4 processes `cd` into your previously installed POET-dir
`<POET_INSTALL_DIR>/bin` and run:
We will continue the above example and start a simulation with *barite_het*,
which simulation files can be found in
`<INSTALL_DIR>/share/poet/barite/barite_het*`. As transport a heterogeneous
diffusion is used. It's a small 2D grid, 2x5 grid, simulating 50 time steps with
a time step size of 100 seconds. To start the simulation with 4 processes `cd`
into your previously installed POET-dir `<POET_INSTALL_DIR>/bin` and run:
```sh
mpirun -n 4 ./poet ../share/poet/bench/dolo/dolo_diffu_inner.R/ output
cp ../share/poet/barite/barite_het* .
mpirun -n 4 ./poet barite_het_rt.R barite_het.rds output
```
After a finished simulation all data generated by POET will be found
@ -200,9 +195,32 @@ produced. This is done by appending the `--dht-snaps=<value>` option. The
resulting call would look like this:
```sh
mpirun -n 4 ./poet --dht --dht-snaps=2 ../share/poet/bench/dolo/dolo_diffu_inner.R/ output
mpirun -n 4 ./poet --dht --dht-snaps=2 barite_het_rt.R barite_het.rds output
```
## Defining a model
In order to provide a model to POET, you need to setup a R script which can then
be used by `poet_init` to generate the simulation input. Which parameters are
required can be found in the
[Wiki](https://git.gfz-potsdam.de/naaice/poet/-/wikis/Initialization). We try to
keep the document up-to-date. However, if you encounter missing information or
need help, please get in touch with us via the issue tracker or E-Mail.
`poet_init` can be used as follows:
```sh
./poet_init [-o, --output output_file] [-s, --setwd] <script.R>
```
where:
- **output** - name of the output file (defaults to the input file name
with the extension `.rds`)
- **setwd** - set the working directory to the directory of the input file (e.g.
to allow relative paths in the input script). However, the output file
will be stored in the directory from which `poet_init` was called.
## About the usage of MPI_Wtime()
Implemented time measurement functions uses `MPI_Wtime()`. Some

56
R_lib/init_r_lib.R Normal file
View File

@ -0,0 +1,56 @@
pqc_to_grid <- function(pqc_in, grid) {
# Convert the input DataFrame to a matrix
pqc_in <- as.matrix(pqc_in)
# Flatten the matrix into a vector
id_vector <- as.numeric(t(grid))
# Find the matching rows in the matrix
row_indices <- match(id_vector, pqc_in[, "ID"])
# Extract the matching rows from pqc_in to size of grid matrix
result_mat <- pqc_in[row_indices, ]
# Convert the result matrix to a data frame
res_df <- as.data.frame(result_mat)
# Remove all columns which only contain NaN
res_df <- res_df[, colSums(is.na(res_df)) != nrow(res_df)]
# Remove row names
rownames(res_df) <- NULL
return(res_df)
}
resolve_pqc_bound <- function(pqc_mat, transport_spec, id) {
df <- as.data.frame(pqc_mat, check.names = FALSE)
value <- df[df$ID == id, transport_spec]
if (is.nan(value)) {
value <- 0
}
return(value)
}
add_missing_transport_species <- function(init_grid, new_names) {
# add 'ID' to new_names front, as it is not a transport species but required
new_names <- c("ID", new_names)
sol_length <- length(new_names)
new_grid <- data.frame(matrix(0, nrow = nrow(init_grid), ncol = sol_length))
names(new_grid) <- new_names
matching_cols <- intersect(names(init_grid), new_names)
# Copy matching columns from init_grid to new_grid
new_grid[, matching_cols] <- init_grid[, matching_cols]
# Add missing columns to new_grid
append_df <- init_grid[, !(names(init_grid) %in% new_names)]
new_grid <- cbind(new_grid, append_df)
return(new_grid)
}

246
R_lib/kin_r_library.R
View File

@ -15,239 +15,105 @@
### this program; if not, write to the Free Software Foundation, Inc., 51
### Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
## Simple function to check file extension. It is needed to check if
## the GridFile is SUM (MUFITS format) or rds/RData
FileExt <- function(x) {
pos <- regexpr("\\.([[:alnum:]]+)$", x)
ifelse(pos > -1L, substring(x, pos + 1L), "")
}
master_init <- function(setup) {
msgm("Process with rank 0 reading GRID properties")
master_init <- function(setup, out_dir, init_field) {
## Setup the directory where we will store the results
verb <- FALSE
if (local_rank == 0) {
verb <- TRUE ## verbosity loading MUFITS results
if (!dir.exists(fileout)) {
dir.create(fileout)
msgm("created directory ", fileout)
} else {
msgm("dir ", fileout, " already exists, I will overwrite!")
}
if (!exists("store_result")) {
msgm("store_result doesn't exist!")
} else {
msgm("store_result is ", store_result)
}
if (!dir.exists(out_dir)) {
dir.create(out_dir)
msgm("created directory ", out_dir)
} else {
msgm("dir ", out_dir, " already exists, I will overwrite!")
}
if (is.null(setup$store_result)) {
msgm("store_result doesn't exist!")
} else {
msgm("store_result is ", setup$store_result)
}
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) {
init_field_out <- paste0(out_dir, "/iter_0.rds")
init_field <- data.frame(init_field, check.names = FALSE)
saveRDS(init_field, file = init_field_out)
msgm("Stored initial field in ", init_field_out)
if (is.null(setup[["out_save"]])) {
setup$out_save <- seq(1, setup$iterations)
}
}
setup$out_dir <- out_dir
return(setup)
}
## This function, called only by master, stores on disk the last
## calculated time step if store_result is TRUE and increments the
## iteration counter
master_iteration_end <- function(setup) {
master_iteration_end <- function(setup, state_T, state_C) {
iter <- setup$iter
# print(iter)
## max digits for iterations
dgts <- as.integer(ceiling(log10(setup$iterations + 1)))
## string format to use in sprintf
dgts <- as.integer(ceiling(log10(setup$maxiter)))
## 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)
)
nameout <- paste0(setup$out_dir, "/iter_", sprintf(fmt = fmt, iter), ".rds")
state_T <- data.frame(state_T, check.names = FALSE)
state_C <- data.frame(state_C, check.names = FALSE)
saveRDS(list(
T = setup$state_T, C = setup$state_C,
simtime = as.integer(setup$simtime),
tr_info = info
T = state_T,
C = state_C,
simtime = as.integer(setup$simulation_time)
), file = nameout)
msgm("results stored in <", nameout, ">")
}
}
msgm("done iteration", iter, "/", setup$maxiter)
## Add last time step to simulation time
setup$simulation_time <- setup$simulation_time + setup$timesteps[iter]
msgm("done iteration", iter, "/", length(setup$timesteps))
setup$iter <- setup$iter + 1
return(setup)
}
## function for the workers to compute chemistry through PHREEQC
slave_chemistry <- function(setup, data) {
base <- setup$base
first <- setup$first
prop <- setup$prop
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])
}
## 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)
}
## 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))
datao <- data[, -indomi]
} 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)) {
red <- cbind(red, rep(data[1, indomi], nrow(red)))
colnames(red)[now + 1] <- names(omit)
ret <- red[, colnames(data)]
} else {
ret <- red
}
rownames(ret) <- NULL
attr(ret, "index") <- inds
return(ret)
}
## 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"))
}
prefix <- paste0("R: ")
cat(paste(prefix, ..., "\n"))
invisible()
}
## Function called by master R process to store on disk all relevant
## parameters for the simulation
StoreSetup <- function(setup) {
StoreSetup <- function(setup, filesim, out_dir) {
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 +148,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
apps/CMakeLists.txt Normal file
View File

@ -0,0 +1,4 @@
file(GLOB INIT_SRCS CONFIGURE_DEPENDS "initializer/*.cpp")
add_executable(poet_initializer ${INIT_SRCS})
target_link_libraries(poet_initializer RRuntime tug)

3
apps/initializer/main.cpp Normal file
View File

@ -0,0 +1,3 @@
#include <Rcpp.h>
int main(int argc, char **argv) {}

43
bench/CMakeLists.txt
View File

@ -1,3 +1,44 @@
function(ADD_BENCH_TARGET TARGET POET_BENCH_LIST RT_FILES OUT_PATH)
set(bench_install_dir share/poet/${OUT_PATH})
# create empty list
set(OUT_FILES_LIST "")
foreach(BENCH_FILE ${${POET_BENCH_LIST}})
get_filename_component(BENCH_NAME ${BENCH_FILE} NAME_WE)
set(OUT_FILE ${CMAKE_CURRENT_BINARY_DIR}/${BENCH_NAME}.rds)
add_custom_command(
OUTPUT ${OUT_FILE}
COMMAND $<TARGET_FILE:poet_init> -o ${OUT_FILE} -s ${CMAKE_CURRENT_SOURCE_DIR}/${BENCH_FILE}
COMMENT "Running poet_init on ${BENCH_FILE}"
DEPENDS poet_init ${CMAKE_CURRENT_SOURCE_DIR}/${BENCH_FILE}
VERBATIM
COMMAND_EXPAND_LISTS
)
list(APPEND OUT_FILES_LIST ${OUT_FILE})
endforeach(BENCH_FILE ${${POET_BENCH_LIST}})
add_custom_target(
${TARGET}
DEPENDS ${OUT_FILES_LIST})
install(FILES ${OUT_FILES_LIST} DESTINATION ${bench_install_dir})
# install all ADD_FILES to the same location
install(FILES ${${RT_FILES}} DESTINATION ${bench_install_dir})
endfunction()
# define target name
set(BENCHTARGET benchmarks)
add_custom_target(${BENCHTARGET} ALL)
add_subdirectory(barite)
add_subdirectory(dolo)
add_subdirectory(surfex)
add_subdirectory(barite)

26
bench/barite/CMakeLists.txt
View File

@ -1,8 +1,20 @@
install(FILES
barite.R
barite_interp_eval.R
barite.pqi
db_barite.dat
DESTINATION
share/poet/bench/barite
# Create a list of files
set(bench_files
barite_200.R
barite_het.R
)
set(runtime_files
barite_200_rt.R
barite_het_rt.R
)
# add_custom_target(barite_bench DEPENDS ${bench_files} ${runtime_files})
ADD_BENCH_TARGET(barite_bench
bench_files
runtime_files
"barite"
)
add_dependencies(${BENCHTARGET} barite_bench)

5
bench/barite/README.org
View File

@ -18,8 +18,9 @@ mpirun -np 4 ./poet --interp barite_interp_eval.R barite_results
* List of Files
- =barite.R=: POET input script for a 20x20 simulation grid
- =barite_interp_eval.R=: POET input script for a 400x200 simulation
- =barite_het.R=: POET input script with homogeneous zones for a 5x2 simulation
grid
- =barite_200.R=: POET input script for a 200x200 simulation
grid
- =db_barite.dat=: PHREEQC database containing the kinetic expressions
for barite and celestite, stripped down from =phreeqc.dat=

149
bench/barite/barite.R
View File

@ -1,149 +0,0 @@
## Time-stamp: "Last modified 2024-01-12 12:39:14 delucia"
database <- normalizePath("../share/poet/bench/barite/db_barite.dat")
input_script <- normalizePath("../share/poet/bench/barite/barite.pqi")
#################################################################
## Section 1 ##
## Grid initialization ##
#################################################################
n <- 20
m <- 20
types <- c("scratch", "phreeqc", "rds")
init_cell <- list(
"H" = 110.0124,
"O" = 55.5087,
"Charge" = -1.216307845207e-09,
"Ba" = 1.E-12,
"Cl" = 2.E-12,
"S(6)" = 6.204727095976e-04,
"Sr" = 6.204727095976e-04,
"Barite" = 0.001,
"Celestite" = 1
)
grid <- list(
n_cells = c(n, m),
s_cells = c(1, 1),
type = types[1],
init_cell = as.data.frame(init_cell, check.names = FALSE),
props = names(init_cell),
database = database,
input_script = input_script
)
##################################################################
## Section 2 ##
## Diffusion parameters and boundary conditions ##
##################################################################
## initial conditions
init_diffu <- list(
"H" = 110.0124,
"O" = 55.5087,
"Charge" = -1.216307845207e-09,
"Ba" = 1.E-12,
"Cl" = 2.E-12,
"S(6)" = 6.204727095976e-04,
"Sr" = 6.204727095976e-04
)
injection_diff <- list(
list(
"H" = 111.0124,
"O" = 55.50622,
"Charge" = -3.336970273297e-08,
"Ba" = 0.1,
"Cl" = 0.2,
"S(6)" = 0,
"Sr" = 0)
)
## diffusion coefficients
alpha_diffu <- c(
"H" = 1E-06,
"O" = 1E-06,
"Charge" = 1E-06,
"Ba" = 1E-06,
"Cl" = 1E-06,
"S(6)" = 1E-06,
"Sr" = 1E-06
)
## vecinj_inner <- list(
## l1 = c(1,20,20),
## l2 = c(2,80,80),
## l3 = c(2,60,80)
## )
boundary <- list(
"N" = c(1,1, rep(0, n-2)),
## "N" = rep(0, n),
"E" = rep(0, n),
"S" = rep(0, n),
# "W" = rep(0, n)
"W" = c(1,1, rep(0, n-2))
)
diffu_list <- names(alpha_diffu)
vecinj <- do.call(rbind.data.frame, injection_diff)
names(vecinj) <- names(init_diffu)
diffusion <- list(
init = as.data.frame(init_diffu, check.names = FALSE),
vecinj = vecinj,
# vecinj_inner = vecinj_inner,
vecinj_index = boundary,
alpha = alpha_diffu
)
#################################################################
## Section 3 ##
## Chemistry module (Phreeqc) ##
#################################################################
## # Needed when using DHT
dht_species <- c(
"H" = 7,
"O" = 7,
"Charge" = 4,
"Ba" = 7,
"Cl" = 7,
"S(6)" = 7,
"Sr" = 7,
"Barite" = 4,
"Celestite" = 4
)
chemistry <- list(
database = database,
input_script = input_script,
dht_species = dht_species
)
#################################################################
## Section 4 ##
## Putting all those things together ##
#################################################################
iterations <- 4
dt <- 100
setup <- list(
grid = grid,
diffusion = diffusion,
chemistry = chemistry,
iterations = iterations,
timesteps = rep(dt, iterations),
store_result = TRUE,
out_save = seq(1, iterations)
)

53
bench/barite/barite.pqi
View File

@ -1,25 +1,32 @@
SELECTED_OUTPUT
-high_precision true
-reset false
-kinetic_reactants Barite Celestite
-saturation_indices Barite Celestite
SOLUTION 1
units mol/kgw
water 1
temperature 25
pH 7
pe 10.799
Ba 0.1
Cl 0.2
S 1e-9
Sr 1e-9
KINETICS 1
Barite
-m 0.001
-parms 50. # reactive surface area
-tol 1e-9
Celestite
-m 1
-parms 10.0 # reactive surface area
-tol 1e-9
units mol/kgw
water 1
temperature 25
pH 7
PURE 1
Celestite 0.0 1
END
RUN_CELLS
-cells 1
COPY solution 1 2
KINETICS 2
Barite
-m 0.001
-parms 50. # reactive surface area
-tol 1e-9
Celestite
-m 1
-parms 10.0 # reactive surface area
-tol 1e-9
END
SOLUTION 3
units mol/kgw
water 1
temperature 25
Ba 0.1
Cl 0.2
END

144
bench/barite/barite_200.R
View File

@ -1,105 +1,39 @@
## Time-stamp: "Last modified 2024-01-12 12:49:03 delucia"
cols <- 200
rows <- 200
database <- normalizePath("../share/poet/bench/barite/db_barite.dat")
input_script <- normalizePath("../share/poet/bench/barite/barite.pqi")
## database <- normalizePath("/home/work/simR/Rphree/poetsims/Sims/Hans/db_barite.dat")
## input_script <- normalizePath("/home/work/simR/Rphree/poetsims/Sims/Hans/barite.pqi")
s_cols <- 1
s_rows <- 1
#################################################################
## Section 1 ##
## Grid initialization ##
#################################################################
grid_def <- matrix(2, nrow = rows, ncol = cols)
n <- 200
m <- 200
init_cell <- list(
"H" = 110.0124,
"O" = 55.5087,
"Charge" = -1.216307845207e-09,
"Ba" = 1.E-12,
"Cl" = 2.E-12,
"S(6)" = 6.204727095976e-04,
"Sr" = 6.204727095976e-04,
"Barite" = 0.001,
"Celestite" = 1
# Define grid configuration for POET model
grid_setup <- list(
pqc_in_file = "./barite.pqi",
pqc_db_file = "./db_barite.dat", # Path to the database file for Phreeqc
grid_def = grid_def, # Definition of the grid, containing IDs according to the Phreeqc input script
grid_size = c(s_rows, s_cols), # Size of the grid in meters
constant_cells = c() # IDs of cells with constant concentration
)
grid <- list(
n_cells = c(n, m),
s_cells = c(1, 1),
type = "scratch",
init_cell = as.data.frame(init_cell, check.names = FALSE),
props = names(init_cell),
database = database,
input_script = input_script
bound_length <- 2
bound_def <- list(
"type" = rep("constant", bound_length),
"sol_id" = rep(3, bound_length),
"cell" = seq(1, bound_length)
)
homogenous_alpha <- 1e-6
##################################################################
## Section 2 ##
## Diffusion parameters and boundary conditions ##
##################################################################
## initial conditions
init_diffu <- list(
"H" = 110.0124,
"O" = 55.5087,
"Charge" = -1.216307845207e-09,
"Ba" = 1.E-12,
"Cl" = 2.E-12,
"S(6)" = 6.204727095976e-04,
"Sr" = 6.204727095976e-04
diffusion_setup <- list(
boundaries = list(
"W" = bound_def,
"N" = bound_def
),
alpha_x = homogenous_alpha,
alpha_y = homogenous_alpha
)
injection_diff <- list(
list(
"H" = 111.0124,
"O" = 55.50622,
"Charge" = -3.336970273297e-08,
"Ba" = 0.1,
"Cl" = 0.2,
"S(6)" = 0,
"Sr" = 0)
)
## diffusion coefficients
alpha_diffu <- c(
"H" = 1E-06,
"O" = 1E-06,
"Charge" = 1E-06,
"Ba" = 1E-06,
"Cl" = 1E-06,
"S(6)" = 1E-06,
"Sr" = 1E-06
)
boundary <- list(
"N" = c(1,1, rep(0, n-2)),
"E" = rep(0, n),
"S" = rep(0, n),
"W" = c(1,1, rep(0, n-2))
)
diffu_list <- names(alpha_diffu)
vecinj <- do.call(rbind.data.frame, injection_diff)
names(vecinj) <- names(init_diffu)
diffusion <- list(
init = as.data.frame(init_diffu, check.names = FALSE),
vecinj = vecinj,
vecinj_index = boundary,
alpha = alpha_diffu
)
#################################################################
## Section 3 ##
## Chemistry module (Phreeqc) ##
#################################################################
## DHT significant digits
dht_species <- c(
"H" = 7,
"O" = 7,
@ -112,27 +46,13 @@ dht_species <- c(
"Celestite" = 4
)
chemistry <- list(
database = database,
input_script = input_script,
dht_species = dht_species
chemistry_setup <- list(
dht_species = dht_species
)
#################################################################
## Section 4 ##
## Putting all those things together ##
#################################################################
iterations <- 50
dt <- 100
# Define a setup list for simulation configuration
setup <- list(
grid = grid,
diffusion = diffusion,
chemistry = chemistry,
iterations = iterations,
timesteps = rep(dt, iterations),
store_result = TRUE,
out_save = seq(1, iterations)
## out_save = c(1, 5, 10, seq(50, iterations, by=50))
Grid = grid_setup, # Parameters related to the grid structure
Diffusion = diffusion_setup, # Parameters related to the diffusion process
Chemistry = chemistry_setup
)

7
bench/barite/barite_200_rt.R Normal file
View File

@ -0,0 +1,7 @@
iterations <- 50
dt <- 100
list(
timesteps = rep(dt, iterations),
store_result = TRUE
)

32
bench/barite/barite_het.R Normal file
View File

@ -0,0 +1,32 @@
grid_def <- matrix(c(2, 3), nrow = 2, ncol = 5)
# Define grid configuration for POET model
grid_setup <- list(
pqc_in_file = "./barite_het.pqi",
pqc_db_file = "./db_barite.dat", # Path to the database file for Phreeqc
grid_def = grid_def, # Definition of the grid, containing IDs according to the Phreeqc input script
grid_size = c(ncol(grid_def), nrow(grid_def)), # Size of the grid in meters
constant_cells = c() # IDs of cells with constant concentration
)
diffusion_setup <- list(
boundaries = list(
"W" = list(
"type" = rep("constant", nrow(grid_def)),
"sol_id" = rep(4, nrow(grid_def)),
"cell" = seq_len(nrow(grid_def))
)
),
alpha_x = 1e-6,
alpha_y = matrix(runif(10, 1e-8, 1e-7),
nrow = nrow(grid_def),
ncol = ncol(grid_def)
)
)
# Define a setup list for simulation configuration
setup <- list(
Grid = grid_setup, # Parameters related to the grid structure
Diffusion = diffusion_setup, # Parameters related to the diffusion process
Chemistry = list()
)

80
bench/barite/barite_het.pqi Normal file
View File

@ -0,0 +1,80 @@
## Initial: everywhere equilibrium with Celestite NB: The aqueous
## solution *resulting* from this calculation is to be used as initial
## state everywhere in the domain
SOLUTION 1
units mol/kgw
water 1
temperature 25
pH 7
pe 4
S(6) 1e-12
Sr 1e-12
Ba 1e-12
Cl 1e-12
PURE 1
Celestite 0.0 1
SAVE SOLUTION 2 # <- phreeqc keyword to store and later reuse these results
END
RUN_CELLS
-cells 1
COPY solution 1 2-3
## Here a 5x2 domain:
|---+---+---+---+---|
-> | 2 | 2 | 2 | 2 | 2 |
4 |---+---+---+---+---|
-> | 3 | 3 | 3 | 3 | 3 |
|---+---+---+---+---|
## East boundary: "injection" of solution 4. North, W, S boundaries: closed
## Here the two distinct zones: nr 2 with kinetics Celestite (initial
## amount is 0, is then allowed to precipitate) and nr 3 with kinetic
## Celestite and Barite (both initially > 0) where the actual
## replacement takes place
#USE SOLUTION 2 <- PHREEQC keyword to reuse the results from previous calculation
KINETICS 2
Celestite
-m 0 # Allowed to precipitate
-parms 10.0
-tol 1e-9
END
#USE SOLUTION 2 <- PHREEQC keyword to reuse the results from previous calculation
KINETICS 3
Barite
-m 0.001
-parms 50.
-tol 1e-9
Celestite
-m 1
-parms 10.0
-tol 1e-9
END
## A BaCl2 solution (nr 4) is "injected" from the left boundary:
SOLUTION 4
units mol/kgw
pH 7
water 1
temp 25
Ba 0.1
Cl 0.2
END
## NB: again, the *result* of the SOLUTION 4 script defines the
## concentration of all elements (+charge, tot H, tot O)
## Ideally, in the initial state SOLUTION 1 we should not have to
## define the 4 elemental concentrations (S(6), Sr, Ba and Cl) but
## obtain them having run once the scripts with the aqueous solution
## resulting from SOLUTION 1 once with KINETICS 2 and once with
## KINETICS 3.
RUN_CELLS
-cells 2-4

4
bench/barite/barite_het_rt.R Normal file
View File

@ -0,0 +1,4 @@
list(
timesteps = rep(50, 100),
store_result = TRUE
)

136
bench/barite/barite_interp_eval.R
View File

@ -1,136 +0,0 @@
## Time-stamp: "Last modified 2024-01-12 11:35:11 delucia"
database <- normalizePath("../share/poet/bench/barite/db_barite.dat")
input_script <- normalizePath("../share/poet/bench/barite/barite.pqi")
#################################################################
## Section 1 ##
## Grid initialization ##
#################################################################
n <- 400
m <- 200
types <- c("scratch", "phreeqc", "rds")
init_cell <- list(
"H" = 110.0124,
"O" = 55.5087,
"Charge" = -1.216307845207e-09,
"Ba" = 1.E-12,
"Cl" = 2.E-12,
"S(6)" = 6.204727095976e-04,
"Sr" = 6.204727095976e-04,
"Barite" = 0.001,
"Celestite" = 1
)
grid <- list(
n_cells = c(n, m),
s_cells = c(1, 1),
type = types[1],
init_cell = as.data.frame(init_cell, check.names = FALSE),
props = names(init_cell),
database = database,
input_script = input_script
)
##################################################################
## Section 2 ##
## Diffusion parameters and boundary conditions ##
##################################################################
## initial conditions
init_diffu <- list(
"H" = 110.0124,
"O" = 55.5087,
"Charge" = -1.216307845207e-09,
"Ba" = 1.E-12,
"Cl" = 2.E-12,
"S(6)" = 6.204727095976e-04,
"Sr" = 6.204727095976e-04
)
injection_diff <- list(
list(
"H" = 111.0124,
"O" = 55.50622,
"Charge" = -3.336970273297e-08,
"Ba" = 0.1,
"Cl" = 0.2,
"S(6)" = 0,
"Sr" = 0)
)
## diffusion coefficients
alpha_diffu <- c(
"H" = 1E-06,
"O" = 1E-06,
"Charge" = 1E-06,
"Ba" = 1E-06,
"Cl" = 1E-06,
"S(6)" = 1E-06,
"Sr" = 1E-06
)
boundary <- list(
"N" = c(1,1, rep(0, n-2)),
"E" = rep(0, n),
"S" = rep(0, n),
"W" = c(1,1, rep(0, n-2))
)
diffu_list <- names(alpha_diffu)
vecinj <- do.call(rbind.data.frame, injection_diff)
names(vecinj) <- names(init_diffu)
diffusion <- list(
init = as.data.frame(init_diffu, check.names = FALSE),
vecinj = vecinj,
vecinj_index = boundary,
alpha = alpha_diffu
)
#################################################################
## Section 3 ##
## Chemistry module (Phreeqc) ##
#################################################################
## # Needed when using DHT
dht_species <- c(
"H" = 7,
"O" = 7,
"Charge" = 4,
"Ba" = 7,
"Cl" = 7,
"S(6)" = 7,
"Sr" = 7,
"Barite" = 4,
"Celestite" = 4
)
chemistry <- list(
database = database,
input_script = input_script,
dht_species = dht_species
)
#################################################################
## Section 4 ##
## Putting all those things together ##
#################################################################
iterations <- 200
dt <- 250
setup <- list(
grid = grid,
diffusion = diffusion,
chemistry = chemistry,
iterations = iterations,
timesteps = rep(dt, iterations),
store_result = TRUE,
out_save = seq(1, iterations)
)

25
bench/dolo/CMakeLists.txt
View File

@ -1,9 +1,18 @@
install(FILES
dolo_diffu_inner.R
dolo_diffu_inner_large.R
dolo_inner.pqi
dolo_interp_long.R
phreeqc_kin.dat
DESTINATION
share/poet/bench/dolo
set(bench_files
dolo_inner_large.R
dolo_interp.R
)
set(runtime_files
dolo_inner_large_rt.R
dolo_interp_rt.R
)
ADD_BENCH_TARGET(
dolo_bench
bench_files
runtime_files
"dolo"
)
add_dependencies(${BENCHTARGET} dolo_bench)

51
bench/dolo/Eval.R
View File

@ -1,51 +0,0 @@
## Time-stamp: "Last modified 2022-12-16 20:26:03 delucia"
source("../../../util/data_evaluation/RFun_Eval.R")
sd <- ReadRTSims("naaice_2d")
sd <- ReadRTSims("Sim2D")
sd <- ReadRTSims("inner")
tim <- readRDS("inner/timings.rds")
simtimes <- sapply(sd, "[","simtime")
## workhorse function to be used with package "animation"
PlotAn <- function(tot, prop, grid, breaks) {
for (step in seq(1, length(tot))) {
snap <- tot[[step]]$C
time <- tot[[step]]$simtime/3600/24
ind <- match(prop, colnames(snap))
Plot2DCellData(snap[,ind], grid=grid, contour=FALSE, breaks=breaks, nlevels=length(breaks), scale=TRUE, main=paste0(prop," after ", time, "days"))
}
}
options(width=110)
library(viridis)
Plot2DCellData(sd$iter_050$C$Cl, nx=1/100, ny=1/100, contour = TRUE,
nlevels = 12, palette = "heat.colors",
rev.palette = TRUE, scale = TRUE, main="Cl")
Plot2DCellData(sd$iter_050$C$Dolomite, nx=100, ny=100, contour = FALSE,
nlevels = 12, palette = "heat.colors",
rev.palette = TRUE, scale = TRUE, )
cairo_pdf("naaice_inner_Dolo.pdf", width=8, height = 6, family="serif")
Plot2DCellData(sd$iter_100$C$Dolomite, nx=100, ny=100, contour = FALSE,
nlevels = 12, palette = "viridis",
rev.palette = TRUE, scale = TRUE, plot.axes = FALSE,
main="2D Diffusion - Dolomite after 2E+4 s (100 iterations)")
dev.off()
cairo_pdf("naaice_inner_Mg.pdf", width=8, height = 6, family="serif")
Plot2DCellData(sd$iter_100$C$Mg, nx=100, ny=100, contour = FALSE,
nlevels = 12, palette = "terrain.colors",
rev.palette = TRUE, scale = TRUE, plot.axes=FALSE,
main="2D Diffusion - Mg after 2E+4 s (100 iterations)")
dev.off()

9
bench/dolo/README.org
View File

@ -18,16 +18,13 @@ mpirun -np 4 ./poet --dht --interp dolo_interp_long.R dolo_interp_long_res
* List of Files
- =dolo_diffu_inner.R=: POET input script for a 100x100 simulation
grid
- =dolo_interp_long.R=: POET input script for a 400x200 simulation
- =dolo_interp.R=: POET input script for a 400x200 simulation
grid
- =dolo_diffu_inner_large.R=: POET input script for a 400x200
simulation grid
- =phreeqc_kin.dat=: PHREEQC database containing the kinetic expressions
for dolomite and celestite, stripped down from =phreeqc.dat=
- =dol.pqi=: PHREEQC input script for the chemical system
# - =dolo.R=: POET input script for a 20x20 simulation grid
# - =dolo_diffu_inner_large.R=: POET input script for a 400x200
# simulation grid
* Chemical system

72
bench/dolo/dol.pqi
View File

@ -1,35 +1,43 @@
SELECTED_OUTPUT
-high_precision true
-reset false
-time
-soln
-temperature true
-water true
-pH
-pe
-totals C Ca Cl Mg
-kinetic_reactants Calcite Dolomite
-equilibrium O2g
SOLUTION 1
units mol/kgw
temp 25.0
water 1
pH 9.91 charge
pe 4.0
C 1.2279E-04
Ca 1.2279E-04
Cl 1E-12
Mg 1E-12
units mol/kgw
water 1
temperature 25
pH 7
pe 4
PURE 1
O2g -0.1675 10
KINETICS 1
Calcite
-m 0.000207
-parms 0.0032
-tol 1e-10
Dolomite
-m 0.0
-parms 0.00032
-tol 1e-10
Calcite 0.0 1
END
RUN_CELLS
-cells 1
COPY solution 1 2
PURE 2
O2g -0.1675 10
KINETICS 2
Calcite
-m 0.000207
-parms 0.05
-tol 1e-10
Dolomite
-m 0.0
-parms 0.005
-tol 1e-10
END
SOLUTION 3
units mol/kgw
water 1
temp 25
Mg 0.001
Cl 0.002
END
SOLUTION 4
units mol/kgw
water 1
temp 25
Mg 0.002
Cl 0.004
END

190
bench/dolo/dolo_diffu_inner.R
View File

@ -1,190 +0,0 @@
## Time-stamp: "Last modified 2023-08-16 17:04:42 mluebke"
database <- normalizePath("../share/poet/bench/dolo/phreeqc_kin.dat")
input_script <- normalizePath("../share/poet/bench/dolo/dolo_inner.pqi")
#################################################################
## Section 1 ##
## Grid initialization ##
#################################################################
n <- 100
m <- 100
types <- c("scratch", "phreeqc", "rds")
init_cell <- list(
"H" = 110.683,
"O" = 55.3413,
"Charge" = -5.0822e-19,
"C(4)" = 1.2279E-4,
"Ca" = 1.2279E-4,
"Cl" = 0,
"Mg" = 0,
"O2g" = 0.499957,
"Calcite" = 2.07e-4,
"Dolomite" = 0
)
grid <- list(
n_cells = c(n, m),
s_cells = c(1, 1),
type = types[1]
)
##################################################################
## Section 2 ##
## Diffusion parameters and boundary conditions ##
##################################################################
## initial conditions
init_diffu <- list(
"H" = 110.683,
"O" = 55.3413,
"Charge" = -5.0822e-19,
"C(4)" = 1.2279E-4,
"Ca" = 1.2279E-4,
"Cl" = 0,
"Mg" = 0
)
## diffusion coefficients
alpha_diffu <- c(
"H" = 1E-6,
"O" = 1E-6,
"Charge" = 1E-6,
"C(4)" = 1E-6,
"Ca" = 1E-6,
"Cl" = 1E-6,
"Mg" = 1E-6
)
## list of boundary conditions/inner nodes
vecinj_diffu <- list(
list(
"H" = 110.683,
"O" = 55.3413,
"Charge" = 1.90431e-16,
"C(4)" = 0,
"Ca" = 0,
"Cl" = 0.002,
"Mg" = 0.001
),
list(
"H" = 110.683,
"O" = 55.3413,
"Charge" = 1.90431e-16,
"C(4)" = 0,
"Ca" = 0.0,
"Cl" = 0.004,
"Mg" = 0.002
)
)
vecinj_inner <- list(
l1 = c(1, 20, 20),
l2 = c(2, 80, 80),
l3 = c(2, 60, 80)
)
boundary <- list(
# "N" = c(1, rep(0, n-1)),
"N" = rep(0, n),
"E" = rep(0, n),
"S" = rep(0, n),
"W" = rep(0, n)
)
diffu_list <- names(alpha_diffu)
vecinj <- do.call(rbind.data.frame, vecinj_diffu)
names(vecinj) <- names(init_diffu)
diffusion <- list(
init = as.data.frame(init_diffu, check.names = FALSE),
vecinj = vecinj,
vecinj_inner = vecinj_inner,
vecinj_index = boundary,
alpha = alpha_diffu
)
#################################################################
## Section 3 ##
## Chemistry module (Phreeqc) ##
#################################################################
## # Needed when using DHT
dht_species <- c(
"H" = 10,
"O" = 10,
"Charge" = 3,
"C(4)" = 5,
"Ca" = 5,
"Cl" = 5,
"Mg" = 5,
"Calcite" = 5,
"Dolomite" = 5
)
check_sign_cal_dol_dht <- function(old, new) {
if ((old["Calcite"] == 0) != (new["Calcite"] == 0)) {
return(TRUE)
}
if ((old["Dolomite"] == 0) != (new["Dolomite"] == 0)) {
return(TRUE)
}
return(FALSE)
}
fuzz_input_dht_keys <- function(input) {
return(input[names(dht_species)])
}
check_sign_cal_dol_interp <- function(to_interp, data_set) {
data_set <- as.data.frame(do.call(rbind, data_set), check.names = FALSE, optional = TRUE)
names(data_set) <- names(dht_species)
cal <- (data_set$Calcite == 0) == (to_interp["Calcite"] == 0)
dol <- (data_set$Dolomite == 0) == (to_interp["Dolomite"] == 0)
cal_dol_same_sig <- cal == dol
return(rev(which(!cal_dol_same_sig)))
}
check_neg_cal_dol <- function(result) {
neg_sign <- (result["Calcite"] <- 0) || (result["Dolomite"] < 0)
return(any(neg_sign))
}
hooks <- list(
dht_fill = check_sign_cal_dol_dht,
dht_fuzz = fuzz_input_dht_keys,
interp_pre_func = check_sign_cal_dol_interp,
interp_post_func = check_neg_cal_dol
)
chemistry <- list(
database = database,
input_script = input_script,
dht_species = dht_species,
hooks = hooks
)
#################################################################
## Section 4 ##
## Putting all those things together ##
#################################################################
iterations <- 10
dt <- 200
setup <- list(
grid = grid,
diffusion = diffusion,
chemistry = chemistry,
iterations = iterations,
timesteps = rep(dt, iterations),
store_result = TRUE
)

190
bench/dolo/dolo_diffu_inner_large.R
View File

@ -1,190 +0,0 @@
## Time-stamp: "Last modified 2023-08-16 17:05:04 mluebke"
database <- normalizePath("../share/poet/bench/dolo/phreeqc_kin.dat")
input_script <- normalizePath("../share/poet/bench/dolo/dolo_inner.pqi")
#################################################################
## Section 1 ##
## Grid initialization ##
#################################################################
n <- 2000
m <- 1000
types <- c("scratch", "phreeqc", "rds")
init_cell <- list(
"H" = 110.683,
"O" = 55.3413,
"Charge" = -5.0822e-19,
"C" = 1.2279E-4,
"Ca" = 1.2279E-4,
"Cl" = 0,
"Mg" = 0,
"O2g" = 0.499957,
"Calcite" = 2.07e-4,
"Dolomite" = 0
)
grid <- list(
n_cells = c(n, m),
s_cells = c(20, 10),
type = types[1]
)
##################################################################
## Section 2 ##
## Diffusion parameters and boundary conditions ##
##################################################################
## initial conditions
init_diffu <- list(
"H" = 110.683,
"O" = 55.3413,
"Charge" = -5.0822e-19,
"C(4)" = 1.2279E-4,
"Ca" = 1.2279E-4,
"Cl" = 0,
"Mg" = 0
)
## diffusion coefficients
alpha_diffu <- c(
"H" = 1E-6,
"O" = 1E-6,
"Charge" = 1E-6,
"C(4)" = 1E-6,
"Ca" = 1E-6,
"Cl" = 1E-6,
"Mg" = 1E-6
)
## list of boundary conditions/inner nodes
vecinj_diffu <- list(
list(
"H" = 110.683,
"O" = 55.3413,
"Charge" = 1.90431e-16,
"C(4)" = 0,
"Ca" = 0,
"Cl" = 0.002,
"Mg" = 0.001
),
list(
"H" = 110.683,
"O" = 55.3413,
"Charge" = 1.90431e-16,
"C(4)" = 0,
"Ca" = 0.0,
"Cl" = 0.004,
"Mg" = 0.002
)
)
vecinj_inner <- list(
l1 = c(1, 400, 200),
l2 = c(2, 1400, 800),
l3 = c(2, 1600, 800)
)
boundary <- list(
# "N" = c(1, rep(0, n-1)),
"N" = rep(0, n),
"E" = rep(0, m),
"S" = rep(0, n),
"W" = rep(0, m)
)
diffu_list <- names(alpha_diffu)
vecinj <- do.call(rbind.data.frame, vecinj_diffu)
names(vecinj) <- names(init_diffu)
diffusion <- list(
init = as.data.frame(init_diffu, check.names = FALSE),
vecinj = vecinj,
vecinj_inner = vecinj_inner,
vecinj_index = boundary,
alpha = alpha_diffu
)
#################################################################
## Section 3 ##
## Chemistry module (Phreeqc) ##
#################################################################
## # Needed when using DHT
dht_species <- c(
"H" = 10,
"O" = 10,
"Charge" = 3,
"C(4)" = 5,
"Ca" = 5,
"Cl" = 5,
"Mg" = 5,
"Calcite" = 5,
"Dolomite" = 5
)
check_sign_cal_dol_dht <- function(old, new) {
if ((old["Calcite"] == 0) != (new["Calcite"] == 0)) {
return(TRUE)
}
if ((old["Dolomite"] == 0) != (new["Dolomite"] == 0)) {
return(TRUE)
}
return(FALSE)
}
fuzz_input_dht_keys <- function(input) {
return(input[names(dht_species)])
}
check_sign_cal_dol_interp <- function(to_interp, data_set) {
data_set <- as.data.frame(do.call(rbind, data_set), check.names = FALSE, optional = TRUE)
names(data_set) <- names(dht_species)
cal <- (data_set$Calcite == 0) == (to_interp["Calcite"] == 0)
dol <- (data_set$Dolomite == 0) == (to_interp["Dolomite"] == 0)
cal_dol_same_sig <- cal == dol
return(rev(which(!cal_dol_same_sig)))
}
check_neg_cal_dol <- function(result) {
neg_sign <- (result["Calcite"] <- 0) || (result["Dolomite"] < 0)
return(any(neg_sign))
}
hooks <- list(
dht_fill = check_sign_cal_dol_dht,
dht_fuzz = fuzz_input_dht_keys,
interp_pre_func = check_sign_cal_dol_interp,
interp_post_func = check_neg_cal_dol
)
chemistry <- list(
database = database,
input_script = input_script,
dht_species = dht_species,
hooks = hooks
)
#################################################################
## Section 4 ##
## Putting all those things together ##
#################################################################
iterations <- 500
dt <- 50
setup <- list(
grid = grid,
diffusion = diffusion,
chemistry = chemistry,
iterations = iterations,
timesteps = rep(dt, iterations),
store_result = TRUE,
out_save = seq(5, iterations, by = 5)
)

28
bench/dolo/dolo_inner.pqi
View File

@ -1,28 +0,0 @@
SELECTED_OUTPUT
-high_precision true
-reset false
-kinetic_reactants Calcite Dolomite
-equilibrium O2g
SOLUTION 1
units mol/kgw
temp 25.0
water 1
pH 9.91 charge
pe 4.0
C 1.2279E-04
Ca 1.2279E-04
Cl 1E-12
Mg 1E-12
PURE 1
O2g -0.1675 10
KINETICS 1
Calcite
-m 0.00020
-parms 0.05
-tol 1e-10
Dolomite
-m 0.0
-parms 0.005
-tol 1e-10
END

BIN
bench/dolo/dolo_inner.rds Normal file
View File

Binary file not shown.

115
bench/dolo/dolo_inner_large.R Normal file
View File

@ -0,0 +1,115 @@
rows <- 2000
cols <- 1000
grid_def <- matrix(2, nrow = rows, ncol = cols)
# Define grid configuration for POET model
grid_setup <- list(
pqc_in_file = "./dol.pqi",
pqc_db_file = "./phreeqc_kin.dat", # Path to the database file for Phreeqc
grid_def = grid_def, # Definition of the grid, containing IDs according to the Phreeqc input script
grid_size = c(cols, rows) / 100, # Size of the grid in meters
constant_cells = c() # IDs of cells with constant concentration
)
bound_size <- 2
diffusion_setup <- list(
inner_boundaries = list(
"row" = c(400, 1400, 1600),
"col" = c(200, 800, 800),
"sol_id" = c(3, 4, 4)
),
alpha_x = 1e-6,
alpha_y = 1e-6
)
check_sign_cal_dol_dht <- function(old, new) {
if ((old["Calcite"] == 0) != (new["Calcite"] == 0)) {
return(TRUE)
}
if ((old["Dolomite"] == 0) != (new["Dolomite"] == 0)) {
return(TRUE)
}
return(FALSE)
}
fuzz_input_dht_keys <- function(input) {
dht_species <- c(
"H" = 3,
"O" = 3,
"Charge" = 3,
"C(4)" = 6,
"Ca" = 6,
"Cl" = 3,
"Mg" = 5,
"Calcite" = 4,
"Dolomite" = 4
)
return(input[names(dht_species)])
}
check_sign_cal_dol_interp <- function(to_interp, data_set) {
dht_species <- c(
"H" = 3,
"O" = 3,
"Charge" = 3,
"C(4)" = 6,
"Ca" = 6,
"Cl" = 3,
"Mg" = 5,
"Calcite" = 4,
"Dolomite" = 4
)
data_set <- as.data.frame(do.call(rbind, data_set), check.names = FALSE, optional = TRUE)
names(data_set) <- names(dht_species)
cal <- (data_set$Calcite == 0) == (to_interp["Calcite"] == 0)
dol <- (data_set$Dolomite == 0) == (to_interp["Dolomite"] == 0)
cal_dol_same_sig <- cal == dol
return(rev(which(!cal_dol_same_sig)))
}
check_neg_cal_dol <- function(result) {
neg_sign <- (result["Calcite"] < 0) || (result["Dolomite"] < 0)
return(neg_sign)
}
# Optional when using Interpolation (example with less key species and custom
# significant digits)
pht_species <- c(
"C(4)" = 3,
"Ca" = 3,
"Mg" = 2,
"Calcite" = 2,
"Dolomite" = 2
)
chemistry_setup <- list(
dht_species = c(
"H" = 3,
"O" = 3,
"Charge" = 3,
"C(4)" = 6,
"Ca" = 6,
"Cl" = 3,
"Mg" = 5,
"Calcite" = 4,
"Dolomite" = 4
),
pht_species = pht_species,
hooks = list(
dht_fill = check_sign_cal_dol_dht,
dht_fuzz = fuzz_input_dht_keys,
interp_pre = check_sign_cal_dol_interp,
interp_post = check_neg_cal_dol
)
)
# Define a setup list for simulation configuration
setup <- list(
Grid = grid_setup, # Parameters related to the grid structure
Diffusion = diffusion_setup, # Parameters related to the diffusion process
Chemistry = chemistry_setup # Parameters related to the chemistry process
)

10
bench/dolo/dolo_inner_large_rt.R Normal file
View File

@ -0,0 +1,10 @@
iterations <- 500
dt <- 50
out_save <- seq(5, iterations, by = 5)
list(
timesteps = rep(dt, iterations),
store_result = TRUE,
out_save = out_save
)

131
bench/dolo/dolo_interp.R Normal file
View File

@ -0,0 +1,131 @@
rows <- 400
cols <- 200
grid_def <- matrix(2, nrow = rows, ncol = cols)
# Define grid configuration for POET model
grid_setup <- list(
pqc_in_file = "./dol.pqi",
pqc_db_file = "./phreeqc_kin.dat", # Path to the database file for Phreeqc
grid_def = grid_def, # Definition of the grid, containing IDs according to the Phreeqc input script
grid_size = c(2.5, 5), # Size of the grid in meters
constant_cells = c() # IDs of cells with constant concentration
)
bound_def_we <- list(
"type" = rep("constant", rows),
"sol_id" = rep(1, rows),
"cell" = seq(1, rows)
)
bound_def_ns <- list(
"type" = rep("constant", cols),
"sol_id" = rep(1, cols),
"cell" = seq(1, cols)
)
diffusion_setup <- list(
boundaries = list(
"W" = bound_def_we,
"E" = bound_def_we,
"N" = bound_def_ns,
"S" = bound_def_ns
),
inner_boundaries = list(
"row" = floor(rows / 2),
"col" = floor(cols / 2),
"sol_id" = c(3)
),
alpha_x = 1e-6,
alpha_y = 1e-6
)
check_sign_cal_dol_dht <- function(old, new) {
if ((old["Calcite"] == 0) != (new["Calcite"] == 0)) {
return(TRUE)
}
if ((old["Dolomite"] == 0) != (new["Dolomite"] == 0)) {
return(TRUE)
}
return(FALSE)
}
fuzz_input_dht_keys <- function(input) {
dht_species <- c(
"H" = 3,
"O" = 3,
"Charge" = 3,
"C(4)" = 6,
"Ca" = 6,
"Cl" = 3,
"Mg" = 5,
"Calcite" = 4,
"Dolomite" = 4
)
return(input[names(dht_species)])
}
check_sign_cal_dol_interp <- function(to_interp, data_set) {
dht_species <- c(
"H" = 3,
"O" = 3,
"Charge" = 3,
"C(4)" = 6,
"Ca" = 6,
"Cl" = 3,
"Mg" = 5,
"Calcite" = 4,
"Dolomite" = 4
)
data_set <- as.data.frame(do.call(rbind, data_set), check.names = FALSE, optional = TRUE)
names(data_set) <- names(dht_species)
cal <- (data_set$Calcite == 0) == (to_interp["Calcite"] == 0)
dol <- (data_set$Dolomite == 0) == (to_interp["Dolomite"] == 0)
cal_dol_same_sig <- cal == dol
return(rev(which(!cal_dol_same_sig)))
}
check_neg_cal_dol <- function(result) {
neg_sign <- (result["Calcite"] < 0) || (result["Dolomite"] < 0)
return(neg_sign)
}
# Optional when using Interpolation (example with less key species and custom
# significant digits)
pht_species <- c(
"C(4)" = 3,
"Ca" = 3,
"Mg" = 2,
"Calcite" = 2,
"Dolomite" = 2
)
chemistry_setup <- list(
dht_species = c(
"H" = 3,
"O" = 3,
"Charge" = 3,
"C(4)" = 6,
"Ca" = 6,
"Cl" = 3,
"Mg" = 5,
"Calcite" = 4,
"Dolomite" = 4
),
pht_species = pht_species,
hooks = list(
dht_fill = check_sign_cal_dol_dht,
dht_fuzz = fuzz_input_dht_keys,
interp_pre = check_sign_cal_dol_interp,
interp_post = check_neg_cal_dol
)
)
# Define a setup list for simulation configuration
setup <- list(
Grid = grid_setup, # Parameters related to the grid structure
Diffusion = diffusion_setup, # Parameters related to the diffusion process
Chemistry = chemistry_setup # Parameters related to the chemistry process
)

204
bench/dolo/dolo_interp_long.R
View File

@ -1,204 +0,0 @@
## Time-stamp: "Last modified 2023-08-16 14:57:25 mluebke"
database <- normalizePath("../share/poet/bench/dolo/phreeqc_kin.dat")
input_script <- normalizePath("../share/poet/bench/dolo/dolo_inner.pqi")
#################################################################
## Section 1 ##
## Grid initialization ##
#################################################################
n <- 400
m <- 200
types <- c("scratch", "phreeqc", "rds")
init_cell <- list(
"H" = 110.683,
"O" = 55.3413,
"Charge" = -5.0822e-19,
"C" = 1.2279E-4,
"Ca" = 1.2279E-4,
"Cl" = 0,
"Mg" = 0,
"O2g" = 0.499957,
"Calcite" = 2.07e-4,
"Dolomite" = 0
)
grid <- list(
n_cells = c(n, m),
s_cells = c(5, 2.5),
type = types[1]
)
##################################################################
## Section 2 ##
## Diffusion parameters and boundary conditions ##
##################################################################
## initial conditions
init_diffu <- list(
"H" = 1.110124E+02,
"O" = 5.550833E+01,
"Charge" = -1.216307659761E-09,
"C(4)" = 1.230067028174E-04,
"Ca" = 1.230067028174E-04,
"Cl" = 0,
"Mg" = 0
)
## diffusion coefficients
alpha_diffu <- c(
"H" = 1E-6,
"O" = 1E-6,
"Charge" = 1E-6,
"C(4)" = 1E-6,
"Ca" = 1E-6,
"Cl" = 1E-6,
"Mg" = 1E-6
)
## list of boundary conditions/inner nodes
vecinj_diffu <- list(
list(
"H" = 1.110124E+02,
"O" = 5.550796E+01,
"Charge" = -3.230390327801E-08,
"C(4)" = 0,
"Ca" = 0,
"Cl" = 0.002,
"Mg" = 0.001
),
list(
"H" = 110.683,
"O" = 55.3413,
"Charge" = 1.90431e-16,
"C(4)" = 0,
"Ca" = 0.0,
"Cl" = 0.004,
"Mg" = 0.002
),
init_diffu
)
vecinj_inner <- list(
l1 = c(1, floor(n / 2), floor(m / 2))
# l2 = c(2,1400,800),
# l3 = c(2,1600,800)
)
boundary <- list(
# "N" = c(1, rep(0, n-1)),
"N" = rep(3, n),
"E" = rep(3, m),
"S" = rep(3, n),
"W" = rep(3, m)
)
diffu_list <- names(alpha_diffu)
vecinj <- do.call(rbind.data.frame, vecinj_diffu)
names(vecinj) <- names(init_diffu)
diffusion <- list(
init = as.data.frame(init_diffu, check.names = FALSE),
vecinj = vecinj,
vecinj_inner = vecinj_inner,
vecinj_index = boundary,
alpha = alpha_diffu
)
#################################################################
## Section 3 ##
## Chemistry module (Phreeqc) ##
#################################################################
## # optional when using DHT
dht_species <- c(
"H" = 3,
"O" = 3,
"Charge" = 3,
"C(4)" = 6,
"Ca" = 6,
"Cl" = 3,
"Mg" = 5,
"Calcite" = 4,
"Dolomite" = 4
)
## # Optional when using Interpolation (example with less key species and custom
## # significant digits)
# pht_species <- c(
# "C(4)" = 3,
# "Ca" = 3,
# "Mg" = 2,
# "Calcite" = 2,
# "Dolomite" = 2
# )
check_sign_cal_dol_dht <- function(old, new) {
if ((old["Calcite"] == 0) != (new["Calcite"] == 0)) {
return(TRUE)
}
if ((old["Dolomite"] == 0) != (new["Dolomite"] == 0)) {
return(TRUE)
}
return(FALSE)
}
fuzz_input_dht_keys <- function(input) {
return(input[names(dht_species)])
}
check_sign_cal_dol_interp <- function(to_interp, data_set) {
data_set <- as.data.frame(do.call(rbind, data_set), check.names = FALSE, optional = TRUE)
names(data_set) <- names(dht_species)
cal <- (data_set$Calcite == 0) == (to_interp["Calcite"] == 0)
dol <- (data_set$Dolomite == 0) == (to_interp["Dolomite"] == 0)
cal_dol_same_sig <- cal == dol
return(rev(which(!cal_dol_same_sig)))
}
check_neg_cal_dol <- function(result) {
neg_sign <- (result["Calcite"] < 0) || (result["Dolomite"] < 0)
return(neg_sign)
}
hooks <- list(
dht_fill = check_sign_cal_dol_dht,
dht_fuzz = fuzz_input_dht_keys,
interp_pre_func = check_sign_cal_dol_interp,
interp_post_func = check_neg_cal_dol
)
chemistry <- list(
database = database,
input_script = input_script,
dht_species = dht_species,
hooks = hooks
# pht_species = pht_species
)
#################################################################
## Section 4 ##
## Putting all those things together ##
#################################################################
iterations <- 20000
dt <- 200
setup <- list(
grid = grid,
diffusion = diffusion,
chemistry = chemistry,
iterations = iterations,
timesteps = rep(dt, iterations),
store_result = TRUE,
out_save = c(1, seq(50, iterations, by = 50))
)

10
bench/dolo/dolo_interp_rt.R Normal file
View File

@ -0,0 +1,10 @@
iterations <- 20000
dt <- 200
out_save <- seq(50, iterations, by = 50)
list(
timesteps = rep(dt, iterations),
store_result = TRUE,
out_save = out_save
)

27
bench/surfex/CMakeLists.txt
View File

@ -1,9 +1,20 @@
install(FILES
ExBase.pqi
ex.R
surfex.R
SurfExBase.pqi
SMILE_2021_11_01_TH.dat
DESTINATION
share/poet/bench/surfex
set(bench_files
# surfex.R
# ex.R
PoetEGU_surfex_500.R
)
set(runtime_files
# surfex_rt.R
# ex_rt.R
PoetEGU_surfex_500_rt.R
)
ADD_BENCH_TARGET(
surfex_bench
bench_files
runtime_files
"surfex"
)
add_dependencies(${BENCHTARGET} surfex_bench)

24
bench/surfex/ExBase.pqi
View File

@ -37,3 +37,27 @@ EXCHANGE 1
Z 0.0012585
Y 0.0009418
END
SOLUTION 2
temp 13
units mol/kgw
C(-4) 2.92438561098248e-21
C(4) 2.65160558871092e-06
Ca 2.89001071336443e-05
Cl 0.000429291158114428
Fe(2) 1.90823391198114e-07
Fe(3) 3.10832423034763e-12
H(0) 2.7888235127385e-15
K 2.5301787e-06
Mg 2.31391999937907e-05
Na 0.00036746969
S(-2) 1.01376078438546e-14
S(2) 1.42247026981542e-19
S(4) 9.49422092568557e-18
S(6) 2.19812504654191e-05
Sr 6.01218519999999e-07
U(4) 4.82255946569383e-12
U(5) 5.49050615347901e-13
U(6) 1.32462838991902e-09
END

40
bench/surfex/PoetEGU_surfex_500.R Normal file
View File

@ -0,0 +1,40 @@
rows <- 500
cols <- 200
grid_left <- matrix(1, nrow = rows, ncol = cols/2)
grid_rght <- matrix(2, nrow = rows, ncol = cols/2)
grid_def <- cbind(grid_left, grid_rght)
# Define grid configuration for POET model
grid_setup <- list(
pqc_in_file = "./SurfexEGU.pqi",
pqc_db_file = "./SMILE_2021_11_01_TH.dat", # Path to the database file for Phreeqc
grid_def = grid_def, # Definition of the grid, containing IDs according to the Phreeqc input script
grid_size = c(10, 4), # Size of the grid in meters
constant_cells = c() # IDs of cells with constant concentration
)
bound_def <- list(
"type" = rep("constant", cols),
"sol_id" = rep(3, cols),
"cell" = seq(1, cols)
)
diffusion_setup <- list(
boundaries = list(
"N" = bound_def
),
alpha_x = matrix(runif(rows*cols))*1e-8,
alpha_y = matrix(runif(rows*cols))*1e-9## ,1e-10
)
chemistry_setup <- list()
# Define a setup list for simulation configuration
setup <- list(
Grid = grid_setup, # Parameters related to the grid structure
Diffusion = diffusion_setup, # Parameters related to the diffusion process
Chemistry = chemistry_setup # Parameters related to the chemistry process
)

11
bench/surfex/PoetEGU_surfex_500_rt.R Normal file
View File

@ -0,0 +1,11 @@
iterations <- 200
dt <- 1000
out_save <- c(1, 2, seq(5, iterations, by=5))
## out_save <- seq(1, iterations)
list(
timesteps = rep(dt, iterations),
store_result = TRUE,
out_save = out_save
)

2
bench/surfex/README.org
View File

@ -20,7 +20,7 @@ mpirun -np 4 ./poet surfex.R surfex_res
- =ex.R=: POET input script for a 100x100 simulation grid, only
exchange
- =ExBase.pqi=: PHREEQC input script for the =ex.R= model
- =surfex.R=: POET input script for a 100x100 simulation grid
- =surfex.R=: POET input script for a 1000x1000 simulation grid
considering both cation exchange and surface complexation
- =SurfExBase.pqi=: PHREEQC input script for the =surfex.R= model
- =SMILE_2021_11_01_TH.dat=: PHREEQC database containing the

24
bench/surfex/SurfExBase.pqi
View File

@ -54,3 +54,27 @@ EXCHANGE 1
Z 0.0012585
Y 0.0009418
END
SOLUTION 2
temp 13
units mol/kgw
C(-4) 2.92438561098248e-21
C(4) 2.65160558871092e-06
Ca 2.89001071336443e-05
Cl 0.000429291158114428
Fe(2) 1.90823391198114e-07
Fe(3) 3.10832423034763e-12
H(0) 2.7888235127385e-15
K 2.5301787e-06
Mg 2.31391999937907e-05
Na 0.00036746969
S(-2) 1.01376078438546e-14
S(2) 1.42247026981542e-19
S(4) 9.49422092568557e-18
S(6) 2.19812504654191e-05
Sr 6.01218519999999e-07
U(4) 4.82255946569383e-12
U(5) 5.49050615347901e-13
U(6) 1.32462838991902e-09
END

108
bench/surfex/SurfexEGU.pqi Normal file
View File

@ -0,0 +1,108 @@
## Time-stamp: "Last modified 2024-04-12 10:59:59 delucia"
## KNOBS
## -logfile false
## -iterations 10000
## -convergence_tolerance 1E-12
## -step_size 2
## -pe_step_size 2
SOLUTION 1 ## Porewater composition Opalinus Clay, WITHOUT radionuclides, AFTER EQUI_PHASES
pe -2.627 ## Eh = -227 mV, Value from Bossart & Thury (2008)-> PC borehole measurement 2003, Eh still decreasing
density 1.01583 ## kg/dm³ = g/cm³
temp 13 ## mean temperature Mont Terri, Bossart & Thury (2008), calculations performed for 25°C
units mol/kgw
## Mean composition
pH 7.064
Na 2.763e-01
Cl 3.228e-01 charge
S(6) 1.653e-02 as SO4
Ca 2.173e-02
Mg 1.740e-02
K 1.902e-03
Sr 4.520e-04
Fe 1.435e-04
U 2.247e-09
SURFACE 1 Opalinus Clay, clay minerals
## calculated with rho_b=2.2903 kg/dm³, poro=0.1662
## 1 dm³ = 13.565641 kg_sed/kg_pw
-equil 1 ## equilibrate with solution 1
-sites_units density ## set unit for binding site density to sites/nm2
-donnan 4.9e-10 ## calculated after Wigger & Van Loon (2018) for ionic strength after equilibration with minerales for pCO2=2.2 log10 bar
# surface density SSA (m2/g) mass (g/kgw)
Kln_aOH 1.155 11. 3798.4 ## Kaolinite 28 wt% (aluminol and silanol sites)
Kln_siOH 1.155
Ill_sOH 0.05 100. 4205.35 ## Illite 31 wt% (weak und strong binding sites)
Ill_wOH 2.26 ## 2 % strong binding sites
Mll_sOH 0.05 100. 813.94 ## Montmorillonite = smektite = 6 wt% (weak und strong binding sites)
Mll_wOH 2.26 ## 2 % strong binding sites
EXCHANGE 1 Exchanger, only illite+montmorillonite
## Illite = 0.225 eq/kg_rock, Montmorillonit = 0.87 eq/kg_rock
-equil 1 ## equilibrate with solution 1
Z 0.9462 ## = Illite
Y 0.70813 ## = Montmorillonite
END
SOLUTION 2 ## Porewater composition Opalinus Clay, WITHOUT radionuclides, AFTER EQUI_PHASES
pe -2.627 ## Eh = -227 mV, Value from Bossart & Thury (2008)-> PC borehole measurement 2003, Eh still decreasing
density 1.01583 ## kg/dm³ = g/cm³
temp 13 ## mean temperature Mont Terri, Bossart & Thury (2008), calculations performed for 25°C
units mol/kgw
## Mean composition
pH 7.064
Na 2.763e-01
Cl 3.228e-01 charge
S(6) 1.653e-02 as SO4
Ca 2.173e-02
Mg 1.740e-02
K 1.902e-03
Sr 4.520e-04
Fe 1.435e-04
U 2.247e-09
SURFACE 2 Opalinus Clay, clay minerals
-equil 2 ## equilibrate with solution 2
-sites_units density ## set unit for binding site density to
## sites/nm2
-donnan 4.9e-10 ## calculated after Wigger & Van Loon (2018)
## for ionic strength after equilibration
## with minerales for pCO2=2.2 log10 bar
## surface density SSA (m2/g) mass (g/kgw)
Kln_aOH 1.155 11. 2798.4 ## Kaolinite 28 wt% (aluminol and silanol sites)
Kln_siOH 1.155
Ill_sOH 0.05 100. 1205.35 ## Illite 31 wt% (weak und strong binding sites)
Ill_wOH 2.26 ## 2 % strong binding sites
Mll_sOH 0.05 100. 113.94 ## Montmorillonite = smektite = 6 wt% (weak und strong binding sites)
Mll_wOH 2.26 ## 2 % strong binding sites
EXCHANGE 2 Exchanger, only illite+montmorillonite
## Illite = 0.225 eq/kg_rock, Montmorillonit = 0.87 eq/kg_rock
-equil 2 ## equilibrate with solution 1
Z 0.5 ## = Illite
Y 0.2 ## = Montmorillonite
END
SOLUTION 3
pe -2.627 ## Eh = -227 mV, Value from Bossart & Thury (2008)-> PC borehole measurement 2003, Eh still decreasing
density 1.01583 ## kg/dm³ = g/cm³
temp 13 ## mean temperature Mont Terri, Bossart & Thury (2008), calculations performed for 25°C
units mol/kgw
## Mean composition
pH 7.064
Na 3.763e-01
Cl 4.228e-01 charge
S(6) 1.653e-02 as SO4
Ca 2.173e-02
Mg 1.740e-02
K 1.902e-03
Sr 4.520e-04
Fe 1.435e-04
U 1e-6
C 1.991e-03
END
RUN_CELLS
END

159
bench/surfex/ex.R
View File

@ -1,140 +1,37 @@
## Time-stamp: "Last modified 2023-08-02 13:59:35 mluebke"
rows <- 100
cols <- 100
database <- normalizePath("./SMILE_2021_11_01_TH.dat")
input_script <- normalizePath("./ExBase.pqi")
grid_def <- matrix(1, nrow = rows, ncol = cols)
cat(paste(":: R This is a test 1\n"))
# Define grid configuration for POET model
grid_setup <- list(
pqc_in_file = "./SurfExBase.pqi",
pqc_db_file = "./SMILE_2021_11_01_TH.dat", # Path to the database file for Phreeqc
grid_def = grid_def, # Definition of the grid, containing IDs according to the Phreeqc input script
grid_size = c(1, 1), # Size of the grid in meters
constant_cells = c() # IDs of cells with constant concentration
)
#################################################################
## Section 1 ##
## Grid initialization ##
#################################################################
bound_def <- list(
"type" = rep("constant", cols),
"sol_id" = rep(2, cols),
"cell" = seq(1, cols)
)
n <- 100
m <- 100
types <- c("scratch", "phreeqc", "rds")
init_cell <- list(H = 1.476571028625e-01,
O = 7.392297218936e-02,
Charge = -1.765225732724e-18,
`C(-4)` = 2.477908970828e-21,
`C(4)` = 2.647623016916e-06,
Ca = 2.889623169138e-05,
Cl = 4.292806181039e-04,
`Fe(2)` =1.908142472666e-07,
`Fe(3)` =3.173306589931e-12,
`H(0)` =2.675642675119e-15,
K = 2.530134809667e-06,
Mg =2.313806319294e-05,
Na =3.674633059628e-04,
`S(-2)` = 8.589766637180e-15,
`S(2)` = 1.205284362720e-19,
`S(4)` = 9.108958772790e-18,
`S(6)` = 2.198092329098e-05,
Sr = 6.012080128154e-07,
`U(4)` = 1.039668623852e-14,
`U(5)` = 1.208394829796e-15,
`U(6)` = 2.976409147150e-12)
grid <- list(
n_cells = c(n, m),
s_cells = c(1, 1),
type = "scratch"
diffusion_setup <- list(
boundaries = list(
"N" = bound_def
),
alpha_x = 1e-6,
alpha_y = 1e-6
)
##################################################################
## Section 2 ##
## Diffusion parameters and boundary conditions ##
##################################################################
vecinj_diffu <- list(
list(H = 0.147659686316291,
O = 0.0739242798146046,
Charge = 7.46361643222701e-20,
`C(-4)` = 2.92438561098248e-21,
`C(4)` = 2.65160558871092e-06,
Ca = 2.89001071336443e-05,
Cl = 0.000429291158114428,
`Fe(2)` = 1.90823391198114e-07,
`Fe(3)` = 3.10832423034763e-12,
`H(0)` = 2.7888235127385e-15,
K = 2.5301787e-06,
Mg = 2.31391999937907e-05,
Na = 0.00036746969,
`S(-2)` = 1.01376078438546e-14,
`S(2)` = 1.42247026981542e-19,
`S(4)` = 9.49422092568557e-18,
`S(6)` = 2.19812504654191e-05,
Sr = 6.01218519999999e-07,
`U(4)` = 4.82255946569383e-12,
`U(5)` = 5.49050615347901e-13,
`U(6)` = 1.32462838991902e-09)
)
vecinj <- do.call(rbind.data.frame, vecinj_diffu)
names(vecinj) <- grid$props
## diffusion coefficients
alpha_diffu <- c(H = 1E-6, O = 1E-6, Charge = 1E-6, `C(-4)` = 1E-6,
`C(4)` = 1E-6, Ca = 1E-6, Cl = 1E-6, `Fe(2)` = 1E-6,
`Fe(3)` = 1E-6, `H(0)` = 1E-6, K = 1E-6, Mg = 1E-6,
Na = 1E-6, `S(-2)` = 1E-6, `S(2)` = 1E-6,
`S(4)` = 1E-6, `S(6)` = 1E-6, Sr = 1E-6,
`U(4)` = 1E-6, `U(5)` = 1E-6, `U(6)` = 1E-6)
## list of boundary conditions/inner nodes
## vecinj_inner <- list(
## list(1,1,1)
## )
boundary <- list(
"N" = rep(1, n),
"E" = rep(0, n),
"S" = rep(0, n),
"W" = rep(0, n)
)
diffu_list <- names(alpha_diffu)
vecinj <- do.call(rbind.data.frame, vecinj_diffu)
names(vecinj) <- names(init_cell)
diffusion <- list(
init = as.data.frame(init_cell, check.names = FALSE),
vecinj = vecinj,
# vecinj_inner = vecinj_inner,
vecinj_index = boundary,
alpha = alpha_diffu
)
#################################################################
## Section 3 ##
## Chemistry module (Phreeqc) ##
#################################################################
chemistry <- list(
database = database,
input_script = input_script
)
#################################################################
## Section 4 ##
## Putting all those things together ##
#################################################################
iterations <- 10
dt <- 200
chemistry_setup <- list()
# Define a setup list for simulation configuration
setup <- list(
grid = grid,
diffusion = diffusion,
chemistry = chemistry,
iterations = iterations,
timesteps = rep(dt, iterations),
store_result = TRUE
)
Grid = grid_setup, # Parameters related to the grid structure
Diffusion = diffusion_setup, # Parameters related to the diffusion process
Chemistry = chemistry_setup # Parameters related to the chemistry process
)

7
bench/surfex/ex_rt.R Normal file
View File

@ -0,0 +1,7 @@
iterations <- 10
dt <- 200
list(
timesteps = rep(dt, iterations),
store_result = TRUE
)

160
bench/surfex/surfex.R
View File

@ -1,141 +1,37 @@
## Time-stamp: "Last modified 2023-08-02 13:59:44 mluebke"
rows <- 1000
cols <- 1000
database <- normalizePath("../share/poet/bench/surfex/SMILE_2021_11_01_TH.dat")
input_script <- normalizePath("../share/poet/bench/surfex/SurfExBase.pqi")
grid_def <- matrix(1, nrow = rows, ncol = cols)
cat(paste(":: R This is a test 1\n"))
# Define grid configuration for POET model
grid_setup <- list(
pqc_in_file = "./SurfExBase.pqi",
pqc_db_file = "./SMILE_2021_11_01_TH.dat", # Path to the database file for Phreeqc
grid_def = grid_def, # Definition of the grid, containing IDs according to the Phreeqc input script
grid_size = c(rows, cols) / 10, # Size of the grid in meters
constant_cells = c() # IDs of cells with constant concentration
)
#################################################################
## Section 1 ##
## Grid initialization ##
#################################################################
bound_def <- list(
"type" = rep("constant", cols),
"sol_id" = rep(2, cols),
"cell" = seq(1, cols)
)
n <- 1000
m <- 1000
types <- c("scratch", "phreeqc", "rds")
init_cell <- list(H = 1.476571028625e-01,
O = 7.392297218936e-02,
Charge = -1.765225732724e-18,
`C(-4)` = 2.477908970828e-21,
`C(4)` = 2.647623016916e-06,
Ca = 2.889623169138e-05,
Cl = 4.292806181039e-04,
`Fe(2)` =1.908142472666e-07,
`Fe(3)` =3.173306589931e-12,
`H(0)` =2.675642675119e-15,
K = 2.530134809667e-06,
Mg =2.313806319294e-05,
Na =3.674633059628e-04,
`S(-2)` = 8.589766637180e-15,
`S(2)` = 1.205284362720e-19,
`S(4)` = 9.108958772790e-18,
`S(6)` = 2.198092329098e-05,
Sr = 6.012080128154e-07,
`U(4)` = 1.039668623852e-14,
`U(5)` = 1.208394829796e-15,
`U(6)` = 2.976409147150e-12)
grid <- list(
n_cells = c(n, m),
s_cells = c(n/10, m/10),
type = "scratch"
diffusion_setup <- list(
boundaries = list(
"N" = bound_def
),
alpha_x = 1e-6,
alpha_y = 1e-6
)
##################################################################
## Section 2 ##
## Diffusion parameters and boundary conditions ##
##################################################################
vecinj_diffu <- list(
list(H = 0.147659686316291,
O = 0.0739242798146046,
Charge = 7.46361643222701e-20,
`C(-4)` = 2.92438561098248e-21,
`C(4)` = 2.65160558871092e-06,
Ca = 2.89001071336443e-05,
Cl = 0.000429291158114428,
`Fe(2)` = 1.90823391198114e-07,
`Fe(3)` = 3.10832423034763e-12,
`H(0)` = 2.7888235127385e-15,
K = 2.5301787e-06,
Mg = 2.31391999937907e-05,
Na = 0.00036746969,
`S(-2)` = 1.01376078438546e-14,
`S(2)` = 1.42247026981542e-19,
`S(4)` = 9.49422092568557e-18,
`S(6)` = 2.19812504654191e-05,
Sr = 6.01218519999999e-07,
`U(4)` = 4.82255946569383e-12,
`U(5)` = 5.49050615347901e-13,
`U(6)` = 1.32462838991902e-09)
)
vecinj <- do.call(rbind.data.frame, vecinj_diffu)
names(vecinj) <- grid$props
## diffusion coefficients
alpha_diffu <- c(H = 1E-6, O = 1E-6, Charge = 1E-6, `C(-4)` = 1E-6,
`C(4)` = 1E-6, Ca = 1E-6, Cl = 1E-6, `Fe(2)` = 1E-6,
`Fe(3)` = 1E-6, `H(0)` = 1E-6, K = 1E-6, Mg = 1E-6,
Na = 1E-6, `S(-2)` = 1E-6, `S(2)` = 1E-6,
`S(4)` = 1E-6, `S(6)` = 1E-6, Sr = 1E-6,
`U(4)` = 1E-6, `U(5)` = 1E-6, `U(6)` = 1E-6)
## list of boundary conditions/inner nodes
## vecinj_inner <- list(
## list(1,1,1)
## )
boundary <- list(
"N" = rep(1, n),
"E" = rep(0, n),
"S" = rep(0, n),
"W" = rep(0, n)
)
diffu_list <- names(alpha_diffu)
vecinj <- do.call(rbind.data.frame, vecinj_diffu)
names(vecinj) <- names(init_cell)
diffusion <- list(
init = as.data.frame(init_cell, check.names = FALSE),
vecinj = vecinj,
# vecinj_inner = vecinj_inner,
vecinj_index = boundary,
alpha = alpha_diffu
)
#################################################################
## Section 3 ##
## Chemistry module (Phreeqc) ##
#################################################################
chemistry <- list(
database = database,
input_script = input_script
)
#################################################################
## Section 4 ##
## Putting all those things together ##
#################################################################
iterations <- 100
dt <- 200
chemistry_setup <- list()
# Define a setup list for simulation configuration
setup <- list(
grid = grid,
diffusion = diffusion,
chemistry = chemistry,
iterations = iterations,
timesteps = rep(dt, iterations),
store_result = TRUE,
out_save = c(5, iterations)
)
Grid = grid_setup, # Parameters related to the grid structure
Diffusion = diffusion_setup, # Parameters related to the diffusion process
Chemistry = chemistry_setup # Parameters related to the chemistry process
)

10
bench/surfex/surfex_rt.R Normal file
View File

@ -0,0 +1,10 @@
iterations <- 100
dt <- 200
out_save <- seq(5, iterations, by = 5)
list(
timesteps = rep(dt, iterations),
store_result = TRUE,
out_save = out_save
)

1
docs/CMakeLists.txt
View File

@ -14,7 +14,6 @@ if(DOXYGEN_FOUND)
doxygen_add_docs(doxygen
${PROJECT_SOURCE_DIR}/src
${PROJECT_SOURCE_DIR}/README.md
${PROJECT_SOURCE_DIR}/docs/Input_Scripts.md
${PROJECT_SOURCE_DIR}/docs/Output.md
COMMENT "Generate html pages")
endif()

86
docs/Input_Scripts.md
View File

@ -1,86 +0,0 @@
# Input Scripts
In the following the expected schemes of the input scripts is described.
Therefore, each section of the input script gets its own chapter. All sections
should return a `list` as results, which are concatenated to one setup list at
the end of the file. All values must have the same name in order to get parsed
by POET.
## Grid initialization
| name | type | description |
|----------------|----------------|-----------------------------------------------------------------------|
| `n_cells` | Numeric Vector | Number of cells in each direction |
| `s_cells` | Numeric Vector | Spatial resolution of grid in each direction |
| `type` | String | Type of initialization, can be set to *scratch*, *phreeqc* or *rds* |
## Diffusion parameters
| name | type | description |
|----------------|----------------------|-------------------------------------------|
| `init` | Named Numeric Vector | Initial state for each diffused species |
| `vecinj` | Data Frame | Defining all boundary conditions row wise |
| `vecinj_inner` | List of Triples | Inner boundaries |
| `vecinj_index` | List of 4 elements | Ghost nodes boundary conditions |
| `alpha` | Named Numeric Vector | Constant alpha for each species |
### Remark on boundary conditions
Each boundary condition should be defined in `vecinj` as a data frame, where one
row holds one boundary condition.
To define inner (constant) boundary conditions, use a list of triples in
`vecinj_inner`, where each triples is defined by $(i,x,y)$. $i$ is defining the
boundary condition, referencing to the row in `vecinj`. $x$ and $y$ coordinates
then defining the position inside the grid.
Ghost nodes are set by `vecinj_index` which is a list containing boundaries for
each celestial direction (**important**: named by `N, E, S, W`). Each direction
is a numeric vector, also representing a row index of the `vecinj` data frame
for each ghost node, starting at the left-most and upper cell respectively. By
setting the boundary condition to $0$, the ghost node is set as closed boundary.
#### Example
Suppose you have a `vecinj` data frame defining 2 boundary conditions and a grid
consisting of $10 \times 10$ grid cells. Grid cell $(1,1)$ should be set to the
first boundary condition and $(5,6)$ to the second. Also, all boundary
conditions for the ghost nodes should be closed. Except the southern boundary,
which should be set to the first boundary condition injection. The following
setup describes how to setup your initial script, where `n` and `m` are the
grids cell count for each direction ($n = m = 10$):
```R
vecinj_inner <- list (
l1 = c(1, 1, 1),
l2 = c(2, 5, 6)
)
vecinj_index <- list(
"N" = rep(0, n),
"E" = rep(0, m),
"S" = rep(1, n),
"W" = rep(0, m)
)
```
## Chemistry parameters
| name | type | description |
|----------------|--------------|----------------------------------------------------------------------------------|
| `database` | String | Path to the Phreeqc database |
| `input_script` | String | Path the the Phreeqc input script |
| `dht_species` | Named Vector | Indicates significant digits to use for each species for DHT rounding. |
| `pht_species` | Named Vector | Indicates significant digits to use for each species for Interpolation rounding. |
## Final setup
| name | type | description |
|----------------|----------------|------------------------------------------------------------|
| `grid` | List | Grid parameter list |
| `diffusion` | List | Diffusion parameter list |
| `chemistry` | List | Chemistry parameter list |
| `iterations` | Numeric Value | Count of iterations |
| `timesteps` | Numeric Vector | $\Delta t$ to use for specific iteration |
| `store_result` | Boolean | Indicates if results should be stored |
| `out_save` | Numeric Vector | *optional:* At which iteration the states should be stored |

60
docs/Output.md
View File

@ -35,34 +35,50 @@ corresponding values can be found in `<OUTPUT_DIRECTORY>/timings.rds`
and possible to read out within a R runtime with
`readRDS("timings.rds")`. There you will find the following values:
| Value | Description |
|--------------------|----------------------------------------------------------------------------|
| simtime | time spent in whole simulation loop without any initialization and cleanup |
| simtime\_transport | measured time in *transport* subroutine |
| simtime\_chemistry | measured time in *chemistry* subroutine (actual parallelized part) |
| Value | Description |
| --------- | -------------------------------------------------------------------------- |
| simtime | time spent in whole simulation loop without any initialization and cleanup |
| chemistry | measured time in *chemistry* subroutine |
| diffusion | measured time in *diffusion* subroutine |
### chemistry subsetting
### Chemistry subsetting
If running parallel there are also measured timings which are subsets of
*simtime\_​chemistry*.
| Value | Description |
| ------------- | --------------------------------------------------------- |
| simtime | overall runtime of chemistry |
| loop | time spent in send/recv loop of master |
| sequential | sequential part of the master (e.g. shuffling field) |
| idle\_master | idling time of the master waiting for workers |
| idle\_worker | idling time (waiting for work from master) of the workers |
| phreeqc\_time | accumulated times for Phreeqc calls of every worker |
| Value | Description |
|-----------------------|-----------------------------------------------------------|
| chemistry\_loop | time spent in send/​recv loop of master |
| chemistry\_sequential | sequential part of master chemistry |
| idle\_master | idling time (waiting for any free worker) of the master |
| idle\_worker | idling time (waiting for work from master) of the workers |
| phreeqc\_time | accumulated times for Phreeqc calls of every worker |
### DHT usage {#DHT-usage}
#### DHT usage
If running in parallel and with activated DHT, two more timings and also
some profiling about the DHT usage are given:
| Value | Description |
|-----------------|---------------------------------------------------------|
| dht\_fill\_time | time to write data to DHT |
| dht\_get\_time | time to retreive data from DHT |
| dh\_hits | count of data points retrieved from DHT |
| dht\_miss | count of misses/count of data points written to DHT |
| --------------- | ------------------------------------------------------- |
| dht\_hits | count of data points retrieved from DHT |
| dht\_evictions | count of data points evicted by another write operation |
| dht\_get\_time | time to retreive data from DHT |
| dht\_fill\_time | time to write data to DHT |
#### Interpolation
If using interpolation, the following values are given:
| Value | Description |
| -------------- | --------------------------------------------------------------------- |
| interp\_w | time spent to write to PHT |
| interp\_r | time spent to read from DHT/PHT/Cache |
| interp\_g | time spent to gather results from DHT |
| interp\_fc | accumulated time spent in interpolation function call |
| interp\_calls | count of interpolations |
| interp\_cached | count of interpolation data sets, which where cached in the local map |
### Diffusion subsetting
| Value | Description |
| --------- | ------------------------------------------ |
| simtime | overall runtime of diffusion |

1
ext/doctest

@ -1 +0,0 @@
Subproject commit ae7a13539fb71f270b87eb2e874fbac80bc8dda2

1
ext/iphreeqc Submodule

@ -0,0 +1 @@
Subproject commit 749fdbc2e9478046bf3f270c70e5800637246712

1
ext/phreeqcrm

@ -1 +0,0 @@
Subproject commit 6ed14c35322a245e3a9776ef262c0ac0eba3b301

2
ext/tug

@ -1 +1 @@
Subproject commit 25855da6b2930559b542bbadb16299932332d6a3
Subproject commit 449647010ab9cdf9e405139f360424a2b21ab3ab

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

17
src/Base/RInsidePOET.hpp
View File

@ -3,6 +3,7 @@
#include <RInside.h>
#include <Rcpp.h>
#include <Rinternals.h>
#include <cstddef>
#include <exception>
#include <optional>
@ -42,6 +43,13 @@ public:
}
}
RHookFunction(SEXP f) {
try {
this->func = Rcpp::Function(f);
} catch (const std::exception &e) {
}
}
template <typename... Args> T operator()(Args... args) const {
if (func.has_value()) {
return (Rcpp::as<T>(this->func.value()(args...)));
@ -50,8 +58,17 @@ public:
}
}
RHookFunction &operator=(const RHookFunction &rhs) {
this->func = rhs.func;
return *this;
}
RHookFunction(const RHookFunction &rhs) { this->func = rhs.func; }
bool isValid() const { return this->func.has_value(); }
SEXP asSEXP() const { return Rcpp::as<SEXP>(this->func.value()); }
private:
std::optional<Rcpp::Function> func;
};

266
src/Base/SimParams.cpp
View File

@ -1,266 +0,0 @@
// 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)
**
** 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 "SimParams.hpp"
#include "../Chemistry/enums.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>
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);
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) {
MSG("Todo");
MSG("See README.md for further information.");
}
return poet::PARSER_HELP;
}
// if positional arguments are missing
else if (!cmdl(2)) {
if (simparams.world_rank == 0) {
ERRMSG("Kin 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 poet::PARSER_ERROR;
}
// collect all parameters which are not known, print them to stderr and return
// with PARSER_ERROR
std::list<std::string> optionsError = validateOptions(cmdl);
if (!optionsError.empty()) {
if (simparams.world_rank == 0) {
ERRMSG("Unrecognized option(s):\n");
for (auto option : optionsError) {
ERRMSG(std::string(option));
}
ERRMSG("Make sure to use available options. Exiting!");
}
return poet::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));
}
}
cmdl(1) >> filesim;
cmdl(2) >> out_dir;
chem_params.dht_outdir = out_dir;
/* distribute information to R runtime */
// if local_rank == 0 then master else worker
R["local_rank"] = simparams.world_rank;
// assign a char* (string) to 'filesim'
R["filesim"] = wrap(filesim);
// 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;
// eval the init string, ignoring any returns
R.parseEvalQ("source(filesim)");
R.parseEvalQ("mysetup <- setup");
this->chem_params.initFromR(R);
return poet::PARSER_OK;
}
void SimParams::initVectorParams(RInside &R) {}
std::list<std::string> SimParams::validateOptions(argh::parser cmdl) {
/* store all unknown parameters here */
std::list<std::string> retList;
/* loop over all flags and compare to given flaglist*/
for (auto &flag : cmdl.flags()) {
if (!(flaglist.find(flag) != flaglist.end()))
retList.push_back(flag);
}
/* and loop also over params and compare to given paramlist */
for (auto &param : cmdl.params()) {
if (!(paramlist.find(param.first) != paramlist.end()))
retList.push_back(param.first);
}
return retList;
}

275
src/Base/SimParams.hpp
View File

@ -1,275 +0,0 @@
/*
** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
** Potsdam)
**
** Copyright (C) 2018-2022 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 PARSER_H
#define PARSER_H
#include "../DataStructures/DataStructures.hpp"
#include "Macros.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>
#include <Rcpp.h>
// BSD-licenced
/** Standard DHT Size. Defaults to 1 GB (1000 MB) */
constexpr uint32_t DHT_SIZE_PER_PROCESS_MB = 1.5E3;
/** Standard work package size */
#define WORK_PACKAGE_SIZE_DEFAULT 32
namespace poet {
enum { PARSER_OK, PARSER_ERROR, PARSER_HELP };
/**
* @brief Defining all simulation parameters
*
*/
typedef struct {
/** Count of processes in MPI_COMM_WORLD */
int world_size;
/** rank of proces in MPI_COMM_WORLD */
int world_rank;
/** indicates if DHT should be used */
bool dht_enabled;
/** apply logarithm to key before rounding */
bool dht_log;
/** indicates if timestep dt differs between iterations */
bool dt_differ;
/** Indicates, when a DHT snapshot should be written */
int dht_snaps;
/** <b>not implemented</b>: How a DHT is distributed over processes */
int dht_strategy;
/** Size of DHt per process in byter */
unsigned int dht_size_per_process;
/** Default significant digit for rounding */
int dht_significant_digits;
/** Default work package size */
unsigned int wp_size;
/** indicates if resulting grid should be stored after every iteration */
bool store_result;
/** indicating whether the progress bar during chemistry simulation should be
* printed or not */
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;
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;
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);
};
/**
* @brief Reads information from program arguments and R runtime
*
* Providing functions to initialize parameters of the simulation using command
* line parameters and parameters from the R runtime. This class will also parse
* arguments from the commandline and decides if argument is known or unknown.
*
* 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);
/**
* @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 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; };
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; };
/**
* @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; };
const auto &getChemParams() const { return chem_params; }
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"};
t_simparams simparams;
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;
std::string filesim;
std::string out_dir;
ChemistryParams chem_params;
};
} // namespace poet
#endif // PARSER_H

96
src/CMakeLists.txt
View File

@ -1,49 +1,85 @@
add_library(poetlib
Base/Grid.cpp
Base/SimParams.cpp
Chemistry/ChemistryModule.cpp
Chemistry/MasterFunctions.cpp
Chemistry/WorkerFunctions.cpp
Chemistry/SurrogateModels/DHT_Wrapper.cpp
Chemistry/SurrogateModels/DHT.c
Chemistry/SurrogateModels/HashFunctions.cpp
Chemistry/SurrogateModels/InterpolationModule.cpp
Chemistry/SurrogateModels/ProximityHashTable.cpp
add_library(POETLib)
target_sources(POETLib
PRIVATE
Init/InitialList.cpp
Init/GridInit.cpp
Init/DiffusionInit.cpp
Init/ChemistryInit.cpp
DataStructures/Field.cpp
Transport/DiffusionModule.cpp
Chemistry/ChemistryModule.cpp
Chemistry/MasterFunctions.cpp
Chemistry/WorkerFunctions.cpp
Chemistry/SurrogateModels/DHT_Wrapper.cpp
Chemistry/SurrogateModels/DHT.c
Chemistry/SurrogateModels/HashFunctions.cpp
Chemistry/SurrogateModels/InterpolationModule.cpp
Chemistry/SurrogateModels/ProximityHashTable.cpp
)
target_link_libraries(poetlib PUBLIC
MPI::MPI_C
${MATH_LIBRARY}
RRuntime
PhreeqcRM
tug
target_include_directories(POETLib PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}")
target_link_libraries(
POETLib
PUBLIC RRuntime
PUBLIC IPhreeqcPOET
PUBLIC tug
PUBLIC MPI::MPI_C
)
target_compile_definitions(poetlib PUBLIC STRICT_R_HEADERS OMPI_SKIP_MPICXX)
# add_library(poetlib
# Base/Grid.cpp
# Base/SimParams.cpp
# Chemistry/ChemistryModule.cpp
# Chemistry/MasterFunctions.cpp
# Chemistry/WorkerFunctions.cpp
# Chemistry/SurrogateModels/DHT_Wrapper.cpp
# Chemistry/SurrogateModels/DHT.c
# Chemistry/SurrogateModels/HashFunctions.cpp
# Chemistry/SurrogateModels/InterpolationModule.cpp
# Chemistry/SurrogateModels/ProximityHashTable.cpp
# Transport/DiffusionModule.cpp
# )
mark_as_advanced(PHREEQCRM_BUILD_MPI PHREEQCRM_DISABLE_OPENMP)
set(PHREEQCRM_DISABLE_OPENMP ON CACHE BOOL "" FORCE)
# target_link_libraries(poetlib PUBLIC
# DataStructures
# Init
# MPI::MPI_C
# ${MATH_LIBRARY}
# RRuntime
# PhreeqcRM
# tug
# )
option(POET_DHT_DEBUG "Build with DHT debug info" OFF)
target_compile_definitions(POETLib PUBLIC STRICT_R_HEADERS OMPI_SKIP_MPICXX)
if(POET_DHT_DEBUG)
target_compile_definitions(poetlib PRIVATE DHT_STATISTICS)
endif()
# mark_as_advanced(PHREEQCRM_BUILD_MPI PHREEQCRM_DISABLE_OPENMP)
# set(PHREEQCRM_DISABLE_OPENMP ON CACHE BOOL "" FORCE)
option(POET_PHT_ADDITIONAL_INFO "Enables additional information in the PHT" OFF)
# option(POET_DHT_DEBUG "Build with DHT debug info" OFF)
if (POET_PHT_ADDITIONAL_INFO)
target_compile_definitions(poetlib PRIVATE POET_PHT_ADD)
endif()
# if(POET_DHT_DEBUG)
# target_compile_definitions(poetlib PRIVATE DHT_STATISTICS)
# endif()
# option(POET_PHT_ADDITIONAL_INFO "Enables additional information in the PHT" OFF)
# if (POET_PHT_ADDITIONAL_INFO)
# target_compile_definitions(poetlib PRIVATE POET_PHT_ADD)
# endif()
file(READ "${PROJECT_SOURCE_DIR}/R_lib/kin_r_library.R" R_KIN_LIB )
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_link_libraries(poet PRIVATE POETLib MPI::MPI_C RRuntime)
target_include_directories(poet PRIVATE "${CMAKE_CURRENT_BINARY_DIR}")
install(TARGETS poet DESTINATION bin)
add_executable(poet_init initializer.cpp)
target_link_libraries(poet_init PRIVATE POETLib RRuntime)
target_include_directories(poet_init PRIVATE "${CMAKE_CURRENT_BINARY_DIR}")
install(TARGETS poet poet_init DESTINATION bin)

23
src/Chemistry/ChemistryDefs.hpp Normal file
View File

@ -0,0 +1,23 @@
#pragma once
#include <cstdint>
#include <vector>
namespace poet {
enum DHT_PROP_TYPES { DHT_TYPE_DEFAULT, DHT_TYPE_CHARGE, DHT_TYPE_TOTAL };
enum CHEMISTRY_OUT_SOURCE { CHEM_PQC, CHEM_DHT, CHEM_INTERP };
struct WorkPackage {
std::size_t size;
std::vector<std::vector<double>> input;
std::vector<std::vector<double>> output;
std::vector<std::uint8_t> mapping;
WorkPackage(std::size_t _size) : size(_size) {
input.resize(size);
output.resize(size);
mapping.resize(size, CHEM_PQC);
}
};
} // namespace poet

239
src/Chemistry/ChemistryModule.cpp
View File

@ -1,24 +1,17 @@
#include "ChemistryModule.hpp"
#include "PhreeqcEngine.hpp"
#include "SurrogateModels/DHT_Wrapper.hpp"
#include "SurrogateModels/Interpolation.hpp"
#include <PhreeqcRM.h>
#include <algorithm>
#include <cassert>
#include <cstddef>
#include <cstdint>
#include <map>
#include <memory>
#include <mpi.h>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
constexpr uint32_t MB_FACTOR = 1E6;
std::vector<double>
inverseDistanceWeighting(const std::vector<std::int32_t> &to_calc,
const std::vector<double> &from,
@ -26,20 +19,21 @@ inverseDistanceWeighting(const std::vector<std::int32_t> &to_calc,
const std::vector<std::vector<double>> &output) {
std::vector<double> results = from;
const std::uint32_t buffer_size = input.size() + 1;
double buffer[buffer_size];
double from_rescaled;
// const std::uint32_t buffer_size = input.size() + 1;
// double buffer[buffer_size];
// double from_rescaled;
const std::uint32_t data_set_n = input.size();
double rescaled[to_calc.size()][data_set_n + 1];
double weights[data_set_n];
// rescaling over all key elements
for (int key_comp_i = 0; key_comp_i < to_calc.size(); key_comp_i++) {
for (std::uint32_t key_comp_i = 0; key_comp_i < to_calc.size();
key_comp_i++) {
const auto output_comp_i = to_calc[key_comp_i];
// rescale input between 0 and 1
for (int point_i = 0; point_i < data_set_n; point_i++) {
for (std::uint32_t point_i = 0; point_i < data_set_n; point_i++) {
rescaled[key_comp_i][point_i] = input[point_i][key_comp_i];
}
@ -50,7 +44,7 @@ inverseDistanceWeighting(const std::vector<std::int32_t> &to_calc,
const double max = *std::max_element(rescaled[key_comp_i],
rescaled[key_comp_i] + data_set_n + 1);
for (int point_i = 0; point_i < data_set_n; point_i++) {
for (std::uint32_t point_i = 0; point_i < data_set_n; point_i++) {
rescaled[key_comp_i][point_i] =
((max - min) != 0
? (rescaled[key_comp_i][point_i] - min) / (max - min)
@ -62,9 +56,10 @@ inverseDistanceWeighting(const std::vector<std::int32_t> &to_calc,
// calculate distances for each data set
double inv_sum = 0;
for (int point_i = 0; point_i < data_set_n; point_i++) {
for (std::uint32_t point_i = 0; point_i < data_set_n; point_i++) {
double distance = 0;
for (int key_comp_i = 0; key_comp_i < to_calc.size(); key_comp_i++) {
for (std::uint32_t key_comp_i = 0; key_comp_i < to_calc.size();
key_comp_i++) {
distance += std::pow(
rescaled[key_comp_i][point_i] - rescaled[key_comp_i][data_set_n], 2);
}
@ -79,7 +74,8 @@ inverseDistanceWeighting(const std::vector<std::int32_t> &to_calc,
// bool has_h = false;
// bool has_o = false;
for (int key_comp_i = 0; key_comp_i < to_calc.size(); key_comp_i++) {
for (std::uint32_t key_comp_i = 0; key_comp_i < to_calc.size();
key_comp_i++) {
const auto output_comp_i = to_calc[key_comp_i];
double key_delta = 0;
@ -91,7 +87,7 @@ inverseDistanceWeighting(const std::vector<std::int32_t> &to_calc,
// has_o = true;
// }
for (int j = 0; j < data_set_n; j++) {
for (std::uint32_t j = 0; j < data_set_n; j++) {
key_delta += weights[j] * output[j][output_comp_i];
}
@ -158,25 +154,30 @@ inverseDistanceWeighting(const std::vector<std::int32_t> &to_calc,
return results;
}
poet::ChemistryModule::ChemistryModule(uint32_t nxyz, uint32_t wp_size,
std::uint32_t maxiter,
const ChemistryParams &chem_param,
MPI_Comm communicator)
: PhreeqcRM(nxyz, 1), group_comm(communicator), wp_size(wp_size),
params(chem_param) {
poet::ChemistryModule::ChemistryModule(
uint32_t wp_size_, const InitialList::ChemistryInit chem_params,
MPI_Comm communicator)
: group_comm(communicator), wp_size(wp_size_), params(chem_params) {
MPI_Comm_rank(communicator, &comm_rank);
MPI_Comm_size(communicator, &comm_size);
MPI_Comm_size(communicator, &this->comm_size);
MPI_Comm_rank(communicator, &this->comm_rank);
this->is_sequential = comm_size == 1;
this->is_master = comm_rank == 0;
this->is_sequential = (this->comm_size == 1);
this->is_master = (this->comm_rank == 0);
this->n_cells = chem_params.total_grid_cells;
if (!is_sequential && is_master) {
MPI_Bcast(&wp_size, 1, MPI_UINT32_T, 0, this->group_comm);
MPI_Bcast(&maxiter, 1, MPI_UINT32_T, 0, this->group_comm);
if (!is_master) {
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_);
// }
}
this->file_pad = std::ceil(std::log10(maxiter + 1));
}
poet::ChemistryModule::~ChemistryModule() {
@ -185,168 +186,6 @@ poet::ChemistryModule::~ChemistryModule() {
}
}
poet::ChemistryModule
poet::ChemistryModule::createWorker(MPI_Comm communicator,
const ChemistryParams &chem_param) {
uint32_t wp_size;
MPI_Bcast(&wp_size, 1, MPI_UINT32_T, 0, communicator);
std::uint32_t maxiter;
MPI_Bcast(&maxiter, 1, MPI_UINT32_T, 0, communicator);
return ChemistryModule(wp_size, wp_size, maxiter, chem_param, communicator);
}
void poet::ChemistryModule::RunInitFile(const std::string &input_script_path) {
if (this->is_master) {
int f_type = CHEM_INIT;
PropagateFunctionType(f_type);
int count = input_script_path.size();
ChemBCast(&count, 1, MPI_INT);
ChemBCast(const_cast<char *>(input_script_path.data()), count, MPI_CHAR);
}
this->RunFile(true, true, false, input_script_path);
this->RunString(true, false, false, "DELETE; -all; PRINT; -warnings 0;");
this->FindComponents();
PhreeqcRM::initializePOET(this->speciesPerModule, this->prop_names);
this->prop_count = prop_names.size();
char exchange = (speciesPerModule[1] == 0 ? -1 : 1);
char kinetics = (speciesPerModule[2] == 0 ? -1 : 1);
char equilibrium = (speciesPerModule[3] == 0 ? -1 : 1);
char surface = (speciesPerModule[4] == 0 ? -1 : 1);
std::vector<int> ic1;
ic1.resize(this->nxyz * 7, -1);
// TODO: hardcoded reaction modules
for (int i = 0; i < nxyz; i++) {
ic1[i] = 1; // Solution 1
ic1[nxyz + i] = equilibrium; // Equilibrium 1
ic1[2 * nxyz + i] = exchange; // Exchange none
ic1[3 * nxyz + i] = surface; // Surface none
ic1[4 * nxyz + i] = -1; // Gas phase none
ic1[5 * nxyz + i] = -1; // Solid solutions none
ic1[6 * nxyz + i] = kinetics; // Kinetics 1
}
this->InitialPhreeqc2Module(ic1);
}
void poet::ChemistryModule::initializeField(const Field &trans_field) {
if (is_master) {
int f_type = CHEM_INIT_SPECIES;
PropagateFunctionType(f_type);
}
std::vector<std::string> essentials_backup{
prop_names.begin() + speciesPerModule[0], prop_names.end()};
std::vector<std::string> new_solution_names{
this->prop_names.begin(), this->prop_names.begin() + speciesPerModule[0]};
if (is_master) {
for (auto &prop : trans_field.GetProps()) {
if (std::find(new_solution_names.begin(), new_solution_names.end(),
prop) == new_solution_names.end()) {
int size = prop.size();
ChemBCast(&size, 1, MPI_INT);
ChemBCast(prop.data(), prop.size(), MPI_CHAR);
new_solution_names.push_back(prop);
}
}
int end = 0;
ChemBCast(&end, 1, MPI_INT);
} else {
constexpr int MAXSIZE = 128;
MPI_Status status;
int recv_size;
char recv_buffer[MAXSIZE];
while (1) {
ChemBCast(&recv_size, 1, MPI_INT);
if (recv_size == 0) {
break;
}
ChemBCast(recv_buffer, recv_size, MPI_CHAR);
recv_buffer[recv_size] = '\0';
new_solution_names.push_back(std::string(recv_buffer));
}
}
// now sort the new values
std::sort(new_solution_names.begin() + 3, new_solution_names.end());
this->SetPOETSolutionNames(new_solution_names);
this->speciesPerModule[0] = new_solution_names.size();
// and append other processes than solutions
std::vector<std::string> new_prop_names = new_solution_names;
new_prop_names.insert(new_prop_names.end(), essentials_backup.begin(),
essentials_backup.end());
std::vector<std::string> old_prop_names{this->prop_names};
this->prop_names = std::move(new_prop_names);
this->prop_count = prop_names.size();
if (is_master) {
this->n_cells = trans_field.GetRequestedVecSize();
std::vector<std::vector<double>> phreeqc_dump(this->nxyz);
this->getDumpedField(phreeqc_dump);
if (is_sequential) {
std::vector<double> init_vec;
for (std::size_t i = 0; i < n_cells; i++) {
init_vec.insert(init_vec.end(), phreeqc_dump[i].begin(),
phreeqc_dump[i].end());
}
const auto tmp_buffer{init_vec};
this->unshuffleField(tmp_buffer, n_cells, prop_count, 1, init_vec);
this->chem_field = Field(n_cells, init_vec, prop_names);
return;
}
std::vector<double> &phreeqc_init = phreeqc_dump[0];
std::vector<std::vector<double>> initial_values;
for (const auto &vec : trans_field.As2DVector()) {
initial_values.push_back(vec);
}
this->base_totals = {initial_values.at(0).at(0),
initial_values.at(1).at(0)};
ChemBCast(base_totals.data(), 2, MPI_DOUBLE);
for (int i = speciesPerModule[0]; i < phreeqc_init.size(); i++) {
std::vector<double> init(n_cells, phreeqc_init[i]);
initial_values.push_back(std::move(init));
}
this->chem_field = Field(n_cells, initial_values, prop_names);
} else {
ChemBCast(base_totals.data(), 2, MPI_DOUBLE);
}
if (this->params.use_dht || this->params.use_interp) {
initializeDHT(this->params.dht_size, this->params.dht_signifs);
setDHTSnapshots(this->params.dht_snaps, this->params.dht_outdir);
setDHTReadFile(this->params.dht_file);
this->dht_enabled = this->params.use_dht;
if (this->params.use_interp) {
initializeInterp(this->params.pht_max_entries, this->params.pht_size,
this->params.interp_min_entries,
this->params.pht_signifs);
this->interp_enabled = this->params.use_interp;
}
}
}
void poet::ChemistryModule::initializeDHT(
uint32_t size_mb, const NamedVector<std::uint32_t> &key_species) {
constexpr uint32_t MB_FACTOR = 1E6;
@ -368,7 +207,7 @@ void poet::ChemistryModule::initializeDHT(
if (key_species.empty()) {
std::vector<std::uint32_t> default_signif(
this->prop_names.size(), DHT_Wrapper::DHT_KEY_SIGNIF_DEFAULT);
this->prop_count, DHT_Wrapper::DHT_KEY_SIGNIF_DEFAULT);
map_copy = NamedVector<std::uint32_t>(this->prop_names, default_signif);
}
@ -382,7 +221,7 @@ void poet::ChemistryModule::initializeDHT(
this->dht = new DHT_Wrapper(dht_comm, dht_size, map_copy, key_indices,
this->prop_names, params.hooks,
this->prop_count, params.use_interp);
this->prop_count, interp_enabled);
this->dht->setBaseTotals(base_totals.at(0), base_totals.at(1));
}
}
@ -393,6 +232,8 @@ inline std::vector<std::int32_t> poet::ChemistryModule::parseDHTSpeciesVec(
std::vector<int32_t> species_indices;
species_indices.reserve(key_species.size());
const auto test = key_species.getNames();
for (const auto &species : key_species.getNames()) {
auto it = std::find(to_compare.begin(), to_compare.end(), species);
if (it == to_compare.end()) {
@ -469,9 +310,9 @@ void poet::ChemistryModule::initializeInterp(
if (key_species.empty()) {
map_copy = this->dht->getKeySpecies();
for (std::size_t i = 0; i < map_copy.size(); i++) {
for (auto i = 0; i < map_copy.size(); i++) {
const std::uint32_t signif =
map_copy[i] - (map_copy[i] > InterpolationModule::COARSE_DIFF
static_cast<std::uint32_t>(map_copy[i]) - (map_copy[i] > InterpolationModule::COARSE_DIFF
? InterpolationModule::COARSE_DIFF
: 0);
map_copy[i] = signif;

120
src/Chemistry/ChemistryModule.hpp
View File

@ -2,16 +2,17 @@
#ifndef CHEMISTRYMODULE_H_
#define CHEMISTRYMODULE_H_
#include "../Base/SimParams.hpp"
#include "../DataStructures/DataStructures.hpp"
#include "DataStructures/Field.hpp"
#include "DataStructures/NamedVector.hpp"
#include "ChemistryDefs.hpp"
#include "Init/InitialList.hpp"
#include "SurrogateModels/DHT_Wrapper.hpp"
#include "SurrogateModels/Interpolation.hpp"
#include <IrmResult.h>
#include <PhreeqcRM.h>
#include "PhreeqcEngine.hpp"
#include <array>
#include <cstddef>
#include <cstdint>
#include <map>
#include <memory>
@ -24,7 +25,7 @@ namespace poet {
* \brief Wrapper around PhreeqcRM to provide POET specific parallelization with
* easy access.
*/
class ChemistryModule : public PhreeqcRM {
class ChemistryModule {
public:
/**
* Creates a new instance of Chemistry module with given grid cell count, work
@ -41,61 +42,67 @@ public:
* \param wp_size Count of grid cells to fill each work package at maximum.
* \param communicator MPI communicator to distribute work in.
*/
ChemistryModule(uint32_t nxyz, uint32_t wp_size, std::uint32_t maxiter,
const ChemistryParams &chem_param, MPI_Comm communicator);
ChemistryModule(uint32_t wp_size,
const InitialList::ChemistryInit chem_params,
MPI_Comm communicator);
/**
* Deconstructor, which frees DHT data structure if used.
*/
~ChemistryModule();
/**
* Parses the input script and extract information needed during runtime.
*
* **Only run by master**.
*
* Database must be loaded beforehand.
*
* \param input_script_path Path to input script to parse.
*/
void RunInitFile(const std::string &input_script_path);
void masterSetField(Field field);
/**
* Run the chemical simulation with parameters set.
*/
void RunCells();
void simulate(double dt);
/**
* Returns the chemical field.
*/
auto GetField() const { return this->chem_field; }
/**
* Returns all known species names, including not only aqueous species, but
* also equilibrium, exchange, surface and kinetic reactants.
*/
auto GetPropNames() const { return this->prop_names; }
// auto GetPropNames() const { return this->prop_names; }
/**
* Return the accumulated runtime in seconds for chemical simulation.
*/
auto GetChemistryTime() const { return this->chem_t; }
/**
* Create a new worker instance inside given MPI communicator.
*
* Wraps communication needed before instanciation can take place.
*
* \param communicator MPI communicator to distribute work in.
*
* \returns A worker instance with fixed work package size.
*/
static ChemistryModule createWorker(MPI_Comm communicator,
const ChemistryParams &chem_param);
void setFilePadding(std::uint32_t maxiter) {
this->file_pad =
static_cast<std::uint8_t>(std::ceil(std::log10(maxiter + 1)));
}
/**
* Default work package size.
*/
static constexpr uint32_t CHEM_DEFAULT_WP_SIZE = 5;
struct SurrogateSetup {
std::vector<std::string> prop_names;
bool dht_enabled;
std::uint32_t dht_size_mb;
bool interp_enabled;
std::uint32_t interp_bucket_size;
std::uint32_t interp_size_mb;
std::uint32_t interp_min_entries;
};
void masterEnableSurrogates(const SurrogateSetup &setup) {
// FIXME: This is a hack to get the prop_names and prop_count from the setup
this->prop_names = setup.prop_names;
this->prop_count = setup.prop_names.size();
this->dht_enabled = setup.dht_enabled;
this->interp_enabled = setup.interp_enabled;
if (this->dht_enabled || this->interp_enabled) {
this->initializeDHT(setup.dht_size_mb, this->params.dht_species);
}
if (this->interp_enabled) {
this->initializeInterp(setup.interp_bucket_size, setup.interp_size_mb,
setup.interp_min_entries,
this->params.interp_species);
}
}
/**
* Intended to alias input parameters for grid initialization with a single
@ -118,17 +125,6 @@ public:
DHT_SNAPS_ITEREND //!< output snapshots after each iteration
};
/**
* **This function has to be run!**
*
* Merge initial values from existing module with the chemistry module and set
* according internal variables.
*
* \param other Field to merge chemistry with. Most likely it is something
* like the diffusion field.
*/
void initializeField(const Field &other);
/**
* **Only called by workers!** Start the worker listening loop.
*/
@ -243,9 +239,7 @@ protected:
const NamedVector<std::uint32_t> &key_species);
enum {
CHEM_INIT,
CHEM_WP_SIZE,
CHEM_INIT_SPECIES,
CHEM_FIELD_INIT,
CHEM_DHT_ENABLE,
CHEM_DHT_SIGNIF_VEC,
CHEM_DHT_SNAPS,
@ -294,7 +288,7 @@ protected:
using worker_list_t = std::vector<struct worker_info_s>;
using workpointer_t = std::vector<double>::iterator;
void MasterRunParallel();
void MasterRunParallel(double dt);
void MasterRunSequential();
void MasterSendPkgs(worker_list_t &w_list, workpointer_t &work_pointer,
@ -320,8 +314,8 @@ protected:
void WorkerPerfToMaster(int type, const struct worker_s &timings);
void WorkerMetricsToMaster(int type);
IRM_RESULT WorkerRunWorkPackage(WorkPackage &work_package, double dSimTime,
double dTimestep);
void WorkerRunWorkPackage(WorkPackage &work_package, double dSimTime,
double dTimestep);
std::vector<uint32_t> CalculateWPSizesVector(uint32_t n_cells,
uint32_t wp_size) const;
@ -344,7 +338,7 @@ protected:
bool is_sequential;
bool is_master;
uint32_t wp_size{CHEM_DEFAULT_WP_SIZE};
uint32_t wp_size;
bool dht_enabled{false};
int dht_snaps_type{DHT_SNAPS_DISABLED};
std::string dht_file_out_dir;
@ -372,7 +366,7 @@ protected:
bool print_progessbar{false};
std::uint32_t file_pad;
std::uint8_t file_pad{1};
double chem_t{0.};
@ -382,11 +376,9 @@ protected:
Field chem_field;
static constexpr int MODULE_COUNT = 5;
const InitialList::ChemistryInit params;
const ChemistryParams &params;
std::array<std::uint32_t, MODULE_COUNT> speciesPerModule{};
std::map<int, std::unique_ptr<PhreeqcEngine>> phreeqc_instances;
};
} // namespace poet

59
src/Chemistry/MasterFunctions.cpp
View File

@ -1,12 +1,9 @@
#include "ChemistryModule.hpp"
#include <PhreeqcRM.h>
#include <algorithm>
#include <cstddef>
#include <cstdint>
#include <mpi.h>
#include <stdexcept>
#include <vector>
std::vector<uint32_t>
@ -308,45 +305,46 @@ inline void poet::ChemistryModule::MasterRecvPkgs(worker_list_t &w_list,
}
}
void poet::ChemistryModule::RunCells() {
void poet::ChemistryModule::simulate(double dt) {
double start_t{MPI_Wtime()};
if (this->is_sequential) {
MasterRunSequential();
return;
}
MasterRunParallel();
MasterRunParallel(dt);
double end_t{MPI_Wtime()};
this->chem_t += end_t - start_t;
}
void poet::ChemistryModule::MasterRunSequential() {
std::vector<double> shuffled_field =
shuffleField(chem_field.AsVector(), n_cells, prop_count, 1);
// std::vector<double> shuffled_field =
// shuffleField(chem_field.AsVector(), n_cells, prop_count, 1);
std::vector<std::vector<double>> input;
for (std::size_t i = 0; i < n_cells; i++) {
input.push_back(
std::vector<double>(shuffled_field.begin() + (i * prop_count),
shuffled_field.begin() + ((i + 1) * prop_count)));
}
// std::vector<std::vector<double>> input;
// for (std::size_t i = 0; i < n_cells; i++) {
// input.push_back(
// std::vector<double>(shuffled_field.begin() + (i * prop_count),
// shuffled_field.begin() + ((i + 1) *
// prop_count)));
// }
this->setDumpedField(input);
PhreeqcRM::RunCells();
this->getDumpedField(input);
// this->setDumpedField(input);
// PhreeqcRM::RunCells();
// this->getDumpedField(input);
shuffled_field.clear();
for (std::size_t i = 0; i < n_cells; i++) {
shuffled_field.insert(shuffled_field.end(), input[i].begin(),
input[i].end());
}
// shuffled_field.clear();
// for (std::size_t i = 0; i < n_cells; i++) {
// shuffled_field.insert(shuffled_field.end(), input[i].begin(),
// input[i].end());
// }
std::vector<double> out_vec{shuffled_field};
unshuffleField(shuffled_field, n_cells, prop_count, 1, out_vec);
chem_field = out_vec;
// std::vector<double> out_vec{shuffled_field};
// unshuffleField(shuffled_field, n_cells, prop_count, 1, out_vec);
// chem_field = out_vec;
}
void poet::ChemistryModule::MasterRunParallel() {
void poet::ChemistryModule::MasterRunParallel(double dt) {
/* declare most of the needed variables here */
double seq_a, seq_b, seq_c, seq_d;
double worker_chemistry_a, worker_chemistry_b;
@ -360,7 +358,6 @@ void poet::ChemistryModule::MasterRunParallel() {
MPI_Barrier(this->group_comm);
double dt = this->PhreeqcRM::GetTimeStep();
static uint32_t iteration = 0;
/* start time measurement of sequential part */
@ -466,3 +463,13 @@ poet::ChemistryModule::CalculateWPSizesVector(uint32_t n_cells,
return wp_sizes_vector;
}
void poet::ChemistryModule::masterSetField(Field field) {
this->chem_field = field;
this->prop_count = field.GetProps().size();
int ftype = CHEM_FIELD_INIT;
PropagateFunctionType(ftype);
ChemBCast(&this->prop_count, 1, MPI_UINT32_T);
}

6
src/Chemistry/SurrogateModels/DHT_Wrapper.cpp
View File

@ -22,13 +22,15 @@
#include "DHT_Wrapper.hpp"
#include "Init/InitialList.hpp"
#include "Rounding.hpp"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <cstring>
#include <iostream>
#include <stdexcept>
#include <vector>
@ -40,7 +42,7 @@ DHT_Wrapper::DHT_Wrapper(MPI_Comm dht_comm, std::uint64_t dht_size,
const NamedVector<std::uint32_t> &key_species,
const std::vector<std::int32_t> &key_indices,
const std::vector<std::string> &_output_names,
const ChemistryParams::Chem_Hook_Functions &_hooks,
const InitialList::ChemistryHookFunctions &_hooks,
uint32_t data_count, bool _with_interp)
: key_count(key_indices.size()), data_count(data_count),
input_key_elements(key_indices), communicator(dht_comm),

17
src/Chemistry/SurrogateModels/DHT_Wrapper.hpp
View File

@ -23,19 +23,18 @@
#ifndef DHT_WRAPPER_H
#define DHT_WRAPPER_H
#include "../../Base/RInsidePOET.hpp"
#include "../../Base/SimParams.hpp"
#include "../../DataStructures/DataStructures.hpp"
#include "../enums.hpp"
#include "HashFunctions.hpp"
#include "Base/RInsidePOET.hpp"
#include "DataStructures/NamedVector.hpp"
#include "Chemistry/ChemistryDefs.hpp"
#include "Init/InitialList.hpp"
#include "LookupKey.hpp"
#include "Rounding.hpp"
#include <array>
#include <cstdint>
#include <limits>
#include <string>
#include <unordered_map>
#include <utility>
#include <vector>
@ -88,7 +87,7 @@ public:
const NamedVector<std::uint32_t> &key_species,
const std::vector<std::int32_t> &key_indices,
const std::vector<std::string> &output_names,
const ChemistryParams::Chem_Hook_Functions &hooks,
const InitialList::ChemistryHookFunctions &hooks,
uint32_t data_count, bool with_interp);
/**
* @brief Destroy the dht wrapper object
@ -260,7 +259,7 @@ private:
const std::vector<std::string> &output_names;
const ChemistryParams::Chem_Hook_Functions &hooks;
const InitialList::ChemistryHookFunctions &hooks;
const bool with_interp;
DHT_ResultObject dht_results;

12
src/Chemistry/SurrogateModels/Interpolation.hpp
View File

@ -3,14 +3,13 @@
#ifndef INTERPOLATION_H_
#define INTERPOLATION_H_
#include "../../Base/SimParams.hpp"
#include "../../DataStructures/DataStructures.hpp"
#include "DataStructures/NamedVector.hpp"
#include "DHT_Wrapper.hpp"
#include "Init/InitialList.hpp"
#include "LookupKey.hpp"
#include "Rounding.hpp"
#include <cassert>
#include <iostream>
#include <list>
#include <memory>
#include <mpi.h>
@ -22,7 +21,6 @@ extern "C" {
}
#include <cstdint>
#include <functional>
#include <unordered_map>
#include <vector>
@ -165,7 +163,7 @@ public:
const NamedVector<std::uint32_t> &interp_key_signifs,
const std::vector<std::int32_t> &dht_key_indices,
const std::vector<std::string> &out_names,
const ChemistryParams::Chem_Hook_Functions &hooks);
const InitialList::ChemistryHookFunctions &hooks);
enum result_status { RES_OK, INSUFFICIENT_DATA, NOT_NEEDED };
@ -261,7 +259,7 @@ private:
return out_key;
}
const ChemistryParams::Chem_Hook_Functions &hooks;
const InitialList::ChemistryHookFunctions &hooks;
const std::vector<std::string> &out_names;
const std::vector<std::string> dht_names;
};

5
src/Chemistry/SurrogateModels/InterpolationModule.cpp
View File

@ -1,8 +1,9 @@
// Time-stamp: "Last modified 2023-08-16 17:02:31 mluebke"
#include "Init/InitialList.hpp"
#include "Interpolation.hpp"
#include "../../DataStructures/DataStructures.hpp"
#include "DHT_Wrapper.hpp"
#include "DataStructures/NamedVector.hpp"
#include "HashFunctions.hpp"
#include "LookupKey.hpp"
#include "Rounding.hpp"
@ -35,7 +36,7 @@ InterpolationModule::InterpolationModule(
const NamedVector<std::uint32_t> &interp_key_signifs,
const std::vector<std::int32_t> &dht_key_indices,
const std::vector<std::string> &_out_names,
const ChemistryParams::Chem_Hook_Functions &_hooks)
const InitialList::ChemistryHookFunctions &_hooks)
: dht_instance(dht), key_signifs(interp_key_signifs),
key_indices(dht_key_indices), min_entries_needed(min_entries_needed),
dht_names(dht.getKeySpecies().getNames()), out_names(_out_names),

77
src/Chemistry/WorkerFunctions.cpp
View File

@ -2,10 +2,9 @@
#include "SurrogateModels/DHT_Wrapper.hpp"
#include "SurrogateModels/Interpolation.hpp"
#include <IrmResult.h>
#include "Chemistry/ChemistryDefs.hpp"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <iomanip>
@ -17,6 +16,7 @@
#include <vector>
namespace poet {
inline std::string get_string(int root, MPI_Comm communicator) {
int count;
MPI_Bcast(&count, 1, MPI_INT, root, communicator);
@ -45,13 +45,8 @@ void poet::ChemistryModule::WorkerLoop() {
PropagateFunctionType(func_type);
switch (func_type) {
case CHEM_INIT: {
RunInitFile(get_string(0, this->group_comm));
break;
}
case CHEM_INIT_SPECIES: {
Field dummy;
initializeField(dummy);
case CHEM_FIELD_INIT: {
ChemBCast(&this->prop_count, 1, MPI_UINT32_T);
break;
}
case CHEM_WORK_LOOP: {
@ -177,9 +172,7 @@ void poet::ChemistryModule::WorkerDoWork(MPI_Status &probe_status,
phreeqc_time_start = MPI_Wtime();
if (WorkerRunWorkPackage(s_curr_wp, current_sim_time, dt) != IRM_OK) {
std::cerr << "Phreeqc error" << std::endl;
};
WorkerRunWorkPackage(s_curr_wp, current_sim_time, dt);
phreeqc_time_end = MPI_Wtime();
@ -285,42 +278,34 @@ void poet::ChemistryModule::WorkerReadDHTDump(
}
}
IRM_RESULT
poet::ChemistryModule::WorkerRunWorkPackage(WorkPackage &work_package,
double dSimTime, double dTimestep) {
// check if we actually need to start phreeqc
std::vector<std::uint32_t> pqc_mapping;
void poet::ChemistryModule::WorkerRunWorkPackage(WorkPackage &work_package,
double dSimTime,
double dTimestep) {
for (std::size_t i = 0; i < work_package.size; i++) {
if (work_package.mapping[i] == CHEM_PQC) {
pqc_mapping.push_back(i);
for (std::size_t wp_id = 0; wp_id < work_package.size; wp_id++) {
if (work_package.mapping[wp_id] != CHEM_PQC) {
continue;
}
auto curr_input = work_package.input[wp_id];
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(std::remove_if(curr_input.begin(), curr_input.end(),
[](double d) { return std::isnan(d); }),
curr_input.end());
phreeqc_instance->runCell(curr_input, dTimestep);
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++];
}
}
}
if (pqc_mapping.empty()) {
return IRM_OK;
}
IRM_RESULT result;
this->PhreeqcRM::setPOETMapping(pqc_mapping);
this->setDumpedField(work_package.input);
this->PhreeqcRM::SetTime(dSimTime);
this->PhreeqcRM::SetTimeStep(dTimestep);
result = this->PhreeqcRM::RunCells();
std::vector<std::vector<double>> output_tmp(work_package.size);
this->getDumpedField(output_tmp);
for (std::size_t i = 0; i < work_package.size; i++) {
if (work_package.mapping[i] == CHEM_PQC) {
work_package.output[i] = output_tmp[i];
}
}
return result;
}
void poet::ChemistryModule::WorkerPerfToMaster(int type,

10
src/Chemistry/enums.hpp
View File

@ -1,10 +0,0 @@
#ifndef ENUMS_H_
#define ENUMS_H_
namespace poet {
enum DHT_PROP_TYPES { DHT_TYPE_DEFAULT, DHT_TYPE_CHARGE, DHT_TYPE_TOTAL };
enum CHEMISTRY_OUT_SOURCE { CHEM_PQC, CHEM_DHT, CHEM_INTERP };
} // namespace poet
#endif // ENUMS_H_

3
src/DataStructures/CMakeLists.txt Normal file
View File

@ -0,0 +1,3 @@
target_sources(POETLib PRIVATE
Field.cpp
)

87
src/DataStructures/Field.cpp
View File

@ -1,6 +1,9 @@
#include "DataStructures.hpp"
#include "Field.hpp"
#include <Rcpp.h>
#include <Rcpp/DataFrame.h>
#include <Rcpp/exceptions.h>
#include <Rcpp/vector/instantiation.h>
#include <Rinternals.h>
#include <cstddef>
#include <cstdint>
@ -117,40 +120,14 @@ poet::FieldColumn &poet::Field::operator[](const std::string &key) {
}
SEXP poet::Field::asSEXP() const {
const std::size_t cols = this->props.size();
Rcpp::List output;
SEXP s_names = PROTECT(Rf_allocVector(STRSXP, cols));
SEXP s_output = PROTECT(Rf_allocVector(VECSXP, cols));
for (std::size_t prop_i = 0; prop_i < this->props.size(); prop_i++) {
const auto &name = this->props[prop_i];
SEXP s_values = PROTECT(Rf_allocVector(REALSXP, this->req_vec_size));
const auto values = this->find(name)->second;
SET_STRING_ELT(s_names, prop_i, Rf_mkChar(name.c_str()));
for (std::size_t i = 0; i < this->req_vec_size; i++) {
REAL(s_values)[i] = values[i];
}
SET_VECTOR_ELT(s_output, prop_i, s_values);
UNPROTECT(1);
for (const auto &elem : this->props) {
const auto map_it = this->find(elem);
output[elem] = Rcpp::wrap(map_it->second);
}
SEXP s_rownames = PROTECT(Rf_allocVector(INTSXP, this->req_vec_size));
for (std::size_t i = 0; i < this->req_vec_size; i++) {
INTEGER(s_rownames)[i] = static_cast<int>(i + 1);
}
Rf_setAttrib(s_output, R_ClassSymbol,
Rf_ScalarString(Rf_mkChar("data.frame")));
Rf_setAttrib(s_output, R_NamesSymbol, s_names);
Rf_setAttrib(s_output, R_RowNamesSymbol, s_rownames);
UNPROTECT(3);
return s_output;
return output;
}
poet::Field &poet::Field::operator=(const FieldColumn &cont_field) {
@ -199,34 +176,26 @@ poet::Field::operator=(const std::vector<FieldColumn> &cont_field) {
void poet::Field::fromSEXP(const SEXP &s_rhs) {
this->clear();
SEXP s_vec = PROTECT(Rf_coerceVector(s_rhs, VECSXP));
SEXP s_names = PROTECT(Rf_getAttrib(s_vec, R_NamesSymbol));
Rcpp::List in_list;
std::size_t cols = static_cast<std::size_t>(Rf_length(s_vec));
this->props.clear();
this->props.reserve(cols);
for (std::size_t i = 0; i < cols; i++) {
const std::string prop_name(CHAR(STRING_ELT(s_names, i)));
this->props.push_back(prop_name);
SEXP s_values = PROTECT(VECTOR_ELT(s_vec, i));
if (i == 0) {
this->req_vec_size = static_cast<std::uint32_t>(Rf_length(s_values));
}
FieldColumn input(this->req_vec_size);
for (std::size_t j = 0; j < this->req_vec_size; j++) {
input[j] = static_cast<double>(REAL(s_values)[j]);
}
UNPROTECT(1);
this->insert({prop_name, input});
try {
in_list = Rcpp::List(s_rhs);
} catch (const Rcpp::exception &e) {
throw std::runtime_error("Input cannot be casted as list.");
}
UNPROTECT(2);
if (in_list.size() == 0) {
return;
}
this->props = Rcpp::as<std::vector<std::string>>(in_list.names());
const Rcpp::NumericVector &in_vec = in_list[0];
this->req_vec_size = static_cast<std::uint32_t>(in_vec.size());
for (const auto &elem : this->props) {
const auto values = Rcpp::as<std::vector<double>>(in_list[elem]);
this->insert({elem, values});
}
}

41
src/DataStructures/DataStructures.hpp → src/DataStructures/Field.hpp
View File

@ -1,54 +1,15 @@
#ifndef DATASTRUCTURES_H_
#define DATASTRUCTURES_H_
#include "../Chemistry/enums.hpp"
#include <Rcpp.h>
#include <cassert>
#include <cinttypes>
#include <cstddef>
#include <iostream>
#include <list>
#include <cstdint>
#include <string>
#include <utility>
#include <vector>
namespace poet {
struct WorkPackage {
std::size_t size;
std::vector<std::vector<double>> input;
std::vector<std::vector<double>> output;
std::vector<std::uint8_t> mapping;
WorkPackage(size_t _size) : size(_size) {
input.resize(size);
output.resize(size);
mapping.resize(size, CHEM_PQC);
}
};
template <typename T> class NamedVector : public Rcpp::NumericVector {
public:
NamedVector() : Rcpp::NumericVector(){};
NamedVector(const std::vector<std::string> &in_names,
const std::vector<T> &in_values)
: Rcpp::NumericVector(Rcpp::wrap(in_values)) {
this->names() = Rcpp::CharacterVector(Rcpp::wrap(in_names));
}
NamedVector(const SEXP &s) : Rcpp::NumericVector(s){};
bool empty() const { return (this->size() == 0); }
std::vector<std::string> getNames() const {
return Rcpp::as<std::vector<std::string>>(this->names());
}
std::vector<T> getValues() const { return Rcpp::as<std::vector<T>>(*this); }
};
using FieldColumn = std::vector<double>;
/**

25
src/DataStructures/NamedVector.hpp Normal file
View File

@ -0,0 +1,25 @@
#pragma once
#include <Rcpp.h>
namespace poet {
template <typename T> class NamedVector : public Rcpp::NumericVector {
public:
NamedVector() : Rcpp::NumericVector(){};
NamedVector(const std::vector<std::string> &in_names,
const std::vector<T> &in_values)
: Rcpp::NumericVector(Rcpp::wrap(in_values)) {
this->names() = Rcpp::CharacterVector(Rcpp::wrap(in_names));
}
NamedVector(const SEXP &s) : Rcpp::NumericVector(s){};
bool empty() const { return (this->size() == 0); }
std::vector<std::string> getNames() const {
return Rcpp::as<std::vector<std::string>>(this->names());
}
std::vector<T> getValues() const { return Rcpp::as<std::vector<T>>(*this); }
};
} // namespace poet

86
src/Init/ChemistryInit.cpp Normal file
View File

@ -0,0 +1,86 @@
#include "InitialList.hpp"
#include <Rcpp.h>
#include <cstddef>
#include <vector>
namespace poet {
void InitialList::initChemistry(const Rcpp::List &chem) {
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"]);
}
if (chem.containsElementNamed("pht_species")) {
this->interp_species = Rcpp::as<NamedVector<uint32_t>>(chem["pht_species"]);
}
if (chem.containsElementNamed("hooks")) {
this->chem_hooks = Rcpp::List(chem["hooks"]);
std::vector<std::string> hook_names = this->chem_hooks.names();
for (const auto &name : hook_names) {
if (this->hook_name_list.find(name) == this->hook_name_list.end()) {
Rcpp::Rcerr << "Unknown chemistry hook: " << name << std::endl;
}
}
}
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 {
ChemistryInit chem_init;
chem_init.total_grid_cells = this->n_cols * this->n_rows;
// chem_init.field_header = this->field_header;
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;
if (this->chem_hooks.size() > 0) {
if (this->chem_hooks.containsElementNamed("dht_fill")) {
chem_init.hooks.dht_fill = this->chem_hooks["dht_fill"];
}
if (this->chem_hooks.containsElementNamed("dht_fuzz")) {
chem_init.hooks.dht_fuzz = this->chem_hooks["dht_fuzz"];
}
if (this->chem_hooks.containsElementNamed("interp_pre")) {
chem_init.hooks.interp_pre = this->chem_hooks["interp_pre"];
}
if (this->chem_hooks.containsElementNamed("interp_post")) {
chem_init.hooks.interp_post = this->chem_hooks["interp_post"];
}
}
return chem_init;
}
} // namespace poet

406
src/Init/DiffusionInit.cpp Normal file
View File

@ -0,0 +1,406 @@
#include <algorithm>
#include <tug/Boundary.hpp>
// leave above Rcpp includes, as eigen seem to have problems with a preceding
// Rcpp include
#include <map>
#include <set>
#include "DataStructures/Field.hpp"
#include "InitialList.hpp"
#include "PhreeqcInit.hpp"
#include <Rcpp/Function.h>
#include <Rcpp/proxy/ProtectedProxy.h>
#include <Rcpp/vector/instantiation.h>
#include <Rinternals.h>
#include <cstddef>
#include <cstdint>
#include <memory>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
namespace poet {
enum SEXP_TYPE { SEXP_IS_LIST, SEXP_IS_VEC };
const std::map<std::uint8_t, std::string> tug_side_mapping = {
{tug::BC_SIDE_RIGHT, "E"},
{tug::BC_SIDE_LEFT, "W"},
{tug::BC_SIDE_TOP, "N"},
{tug::BC_SIDE_BOTTOM, "S"}};
static std::vector<TugType> colMajToRowMaj(const Rcpp::NumericVector &vec,
std::uint32_t n_cols,
std::uint32_t n_rows) {
if (vec.size() == 1) {
return std::vector<TugType>(n_cols * n_rows, vec[0]);
} else {
if (vec.size() != n_cols * n_rows) {
throw std::runtime_error("Alpha matrix does not match grid dimensions");
}
std::vector<TugType> alpha(n_cols * n_rows);
for (std::uint32_t i = 0; i < n_cols; i++) {
for (std::uint32_t j = 0; j < n_rows; j++) {
alpha[i * n_rows + j] = vec[j * n_cols + i];
}
}
return alpha;
}
}
static std::vector<std::string>
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) {
std::vector<std::string> transport_names = old_trans_names;
std::set<int> constant_pqc_ids;
for (const auto &side : tug_side_mapping) {
if (!boundaries_list.containsElementNamed(side.second.c_str())) {
continue;
}
Rcpp::List mapping = boundaries_list[side.second];
const Rcpp::NumericVector cells = mapping["cell"];
const Rcpp::NumericVector values = mapping["sol_id"];
const Rcpp::CharacterVector type_str = mapping["type"];
if (cells.size() != values.size()) {
throw std::runtime_error("Boundary [" + side.second +
"] cells and values are not the same "
"length");
}
for (auto i = 0; i < cells.size(); i++) {
if (type_str[i] == "constant") {
constant_pqc_ids.insert(static_cast<std::size_t>(values[i]));
}
}
}
if (inner_boundaries.size() > 0) {
const Rcpp::NumericVector values = inner_boundaries["sol_id"];
for (auto i = 0; i < values.size(); i++) {
constant_pqc_ids.insert(static_cast<std::size_t>(values[i]));
}
}
if (!constant_pqc_ids.empty()) {
constexpr std::size_t keep_h_o_charge = 3;
for (const auto &pqc_id : constant_pqc_ids) {
const auto solution_names = phreeqc->getSolutionNames(pqc_id);
// add those strings which are not already in the transport_names
for (const auto &name : solution_names) {
if (std::find(transport_names.begin(), transport_names.end(), name) ==
transport_names.end()) {
auto position =
std::lower_bound(transport_names.begin() + keep_h_o_charge,
transport_names.end(), name);
transport_names.insert(position, name);
}
}
}
}
return transport_names;
}
static Rcpp::List
extend_initial_grid(const Rcpp::List &initial_grid,
const std::vector<std::string> &transport_names) {
// std::vector<std::string> names_to_add(transport_names.begin() + old_size,
// transport_names.end());
Rcpp::Function extend_grid_R("add_missing_transport_species");
return extend_grid_R(initial_grid, Rcpp::wrap(transport_names));
}
std::pair<Rcpp::List, Rcpp::List>
InitialList::resolveBoundaries(const Rcpp::List &boundaries_list,
const Rcpp::List &inner_boundaries) {
Rcpp::List bound_list;
Rcpp::List inner_bound;
Rcpp::Function resolve_R("resolve_pqc_bound");
const std::size_t old_transport_size = this->transport_names.size();
this->transport_names = extend_transport_names(
this->phreeqc, boundaries_list, inner_boundaries, this->transport_names);
const std::size_t new_transport_size = this->transport_names.size();
if (old_transport_size != new_transport_size) {
this->initial_grid =
extend_initial_grid(this->initial_grid, this->transport_names);
}
for (const auto &species : this->transport_names) {
Rcpp::List spec_list;
for (const auto &side : tug_side_mapping) {
Rcpp::NumericVector c_type =
(side.first == tug::BC_SIDE_RIGHT || side.first == tug::BC_SIDE_LEFT)
? Rcpp::NumericVector(n_rows, tug::BC_TYPE_CLOSED)
: Rcpp::NumericVector(n_cols, tug::BC_TYPE_CLOSED);
Rcpp::NumericVector c_value =
(side.first == tug::BC_SIDE_RIGHT || side.first == tug::BC_SIDE_LEFT)
? Rcpp::NumericVector(n_rows, 0)
: Rcpp::NumericVector(n_cols, 0);
if (boundaries_list.containsElementNamed(side.second.c_str())) {
const Rcpp::List mapping = boundaries_list[side.second];
const Rcpp::NumericVector cells = mapping["cell"];
const Rcpp::NumericVector values = mapping["sol_id"];
const Rcpp::CharacterVector type_str = mapping["type"];
if (cells.size() != values.size()) {
throw std::runtime_error(
"Boundary [" + side.second +
"] cells and values are not the same length");
}
for (auto i = 0; i < cells.size(); i++) {
const auto c_id = cells[i] - 1;
if (type_str[i] == "closed") {
c_type[c_id] = tug::BC_TYPE_CLOSED;
} else if (type_str[i] == "constant") {
c_type[c_id] = tug::BC_TYPE_CONSTANT;
c_value[c_id] = Rcpp::as<TugType>(
resolve_R(this->phreeqc_mat, Rcpp::wrap(species), values[i]));
} else {
throw std::runtime_error("Unknown boundary type");
}
}
}
spec_list[side.second] = Rcpp::List::create(
Rcpp::Named("type") = c_type, Rcpp::Named("value") = c_value);
}
bound_list[species] = spec_list;
if (inner_boundaries.size() > 0) {
const Rcpp::NumericVector &inner_row_vec = inner_boundaries["row"];
const Rcpp::NumericVector &inner_col_vec = inner_boundaries["col"];
const Rcpp::NumericVector &inner_pqc_id_vec = inner_boundaries["sol_id"];
std::vector<std::uint32_t> rows;
std::vector<std::uint32_t> cols;
std::vector<TugType> c_value;
if (inner_row_vec.size() != inner_col_vec.size() ||
inner_row_vec.size() != inner_pqc_id_vec.size()) {
throw std::runtime_error(
"Inner boundary vectors are not the same length");
}
for (std::size_t i = 0; i < inner_row_vec.size(); i++) {
rows.push_back(inner_row_vec[i] - 1);
cols.push_back(inner_col_vec[i] - 1);
c_value.push_back(Rcpp::as<TugType>(resolve_R(
this->phreeqc_mat, Rcpp::wrap(species), inner_pqc_id_vec[i])));
}
inner_bound[species] =
Rcpp::List::create(Rcpp::Named("row") = Rcpp::wrap(rows),
Rcpp::Named("col") = Rcpp::wrap(cols),
Rcpp::Named("value") = Rcpp::wrap(c_value));
}
}
return std::make_pair(bound_list, inner_bound);
}
static inline SEXP_TYPE get_datatype(const SEXP &input) {
Rcpp::Function check_list("is.list");
if (Rcpp::as<bool>(check_list(input))) {
return SEXP_IS_LIST;
} else {
return SEXP_IS_VEC;
}
}
static Rcpp::List parseAlphas(const SEXP &input,
const std::vector<std::string> &transport_names,
std::uint32_t n_cols, std::uint32_t n_rows) {
Rcpp::List out_list;
SEXP_TYPE input_type = get_datatype(input);
switch (input_type) {
case SEXP_IS_LIST: {
const Rcpp::List input_list(input);
if (input_list.size() != transport_names.size()) {
throw std::runtime_error(
"Length of alphas list does not match number of transport names");
}
for (const auto &name : transport_names) {
if (!input_list.containsElementNamed(name.c_str())) {
throw std::runtime_error("Alphas list does not contain transport name");
}
const Rcpp::NumericVector &alpha_col_order_vec = input_list[name];
out_list[name] =
Rcpp::wrap(colMajToRowMaj(alpha_col_order_vec, n_cols, n_rows));
}
break;
}
case SEXP_IS_VEC: {
const Rcpp::NumericVector alpha(input);
for (const auto &name : transport_names) {
out_list[name] = Rcpp::wrap(colMajToRowMaj(alpha, n_cols, n_rows));
}
break;
}
default:
throw std::runtime_error("Unknown data type");
}
return out_list;
}
void InitialList::initDiffusion(const Rcpp::List &diffusion_input) {
Rcpp::List boundaries;
Rcpp::List inner_boundaries;
if (diffusion_input.containsElementNamed(
DIFFU_MEMBER_STR(DiffusionMembers::BOUNDARIES))) {
boundaries =
diffusion_input[DIFFU_MEMBER_STR(DiffusionMembers::BOUNDARIES)];
}
if (diffusion_input.containsElementNamed(
DIFFU_MEMBER_STR(DiffusionMembers::INNER_BOUNDARIES))) {
inner_boundaries =
diffusion_input[DIFFU_MEMBER_STR(DiffusionMembers::INNER_BOUNDARIES)];
}
const SEXP &alpha_x =
diffusion_input[DIFFU_MEMBER_STR(DiffusionMembers::ALPHA_X)];
const SEXP &alpha_y =
diffusion_input[DIFFU_MEMBER_STR(DiffusionMembers::ALPHA_Y)];
const auto resolved_boundaries =
resolveBoundaries(boundaries, inner_boundaries);
this->boundaries = resolved_boundaries.first;
this->inner_boundaries = resolved_boundaries.second;
this->alpha_x =
parseAlphas(alpha_x, this->transport_names, this->n_cols, this->n_rows);
this->alpha_y =
parseAlphas(alpha_y, this->transport_names, this->n_cols, this->n_rows);
}
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 (static_cast<tug::BC_TYPE>(type[i]) == tug::BC_TYPE_CONSTANT) {
bc.setBoundaryElementConstant(static_cast<tug::BC_SIDE>(tug_index), i,
value[i]);
}
}
}
map[name] = bc.serialize();
}
return map;
}
static InitialList::DiffusionInit::InnerBoundaryMap
RcppListToInnerBoundaryMap(const std::vector<std::string> &trans_names,
const Rcpp::List &inner_bound_list,
std::uint32_t n_cols, std::uint32_t n_rows) {
InitialList::DiffusionInit::InnerBoundaryMap map;
if (inner_bound_list.size() == 0) {
return map;
}
for (const auto &name : trans_names) {
const Rcpp::List &conc_list = inner_bound_list[name];
std::map<std::pair<std::uint32_t, std::uint32_t>, TugType> inner_bc;
const Rcpp::NumericVector &row = conc_list["row"];
const Rcpp::NumericVector &col = conc_list["col"];
const Rcpp::NumericVector &value = conc_list["value"];
if (row.size() != col.size() || row.size() != value.size()) {
throw std::runtime_error(
"Inner boundary vectors are not the same length");
}
for (std::size_t i = 0; i < row.size(); i++) {
inner_bc[std::make_pair(row[i], col[i])] = value[i];
}
map[name] = inner_bc;
}
return map;
}
InitialList::DiffusionInit InitialList::getDiffusionInit() const {
DiffusionInit diff_init;
diff_init.dim = this->dim;
diff_init.n_cols = this->n_cols;
diff_init.n_rows = this->n_rows;
diff_init.s_cols = this->s_cols;
diff_init.s_rows = this->s_rows;
diff_init.constant_cells = this->constant_cells;
diff_init.transport_names = this->transport_names;
diff_init.boundaries = RcppListToBoundaryMap(
this->transport_names, this->boundaries, this->n_cols, this->n_rows);
diff_init.inner_boundaries =
RcppListToInnerBoundaryMap(this->transport_names, this->inner_boundaries,
this->n_cols, this->n_rows);
diff_init.alpha_x = Field(this->alpha_x);
diff_init.alpha_y = Field(this->alpha_y);
return diff_init;
}
} // namespace poet

190
src/Init/GridInit.cpp Normal file
View File

@ -0,0 +1,190 @@
#include "InitialList.hpp"
#include "PhreeqcInit.hpp"
#include <RInside.h>
#include <Rcpp/Function.h>
#include <Rcpp/vector/instantiation.h>
#include <cstdint>
#include <map>
#include <memory>
#include <regex>
#include <sstream>
#include <string>
#include <vector>
namespace poet {
static Rcpp::NumericMatrix
pqcScriptToGrid(std::unique_ptr<PhreeqcInit> &phreeqc, RInside &R) {
PhreeqcInit::PhreeqcMat phreeqc_mat = phreeqc->getPhreeqcMat();
// add "id" to the front of the column names
const std::size_t ncols = phreeqc_mat.names.size();
const std::size_t nrows = phreeqc_mat.values.size();
phreeqc_mat.names.insert(phreeqc_mat.names.begin(), "ID");
Rcpp::NumericMatrix mat(nrows, ncols + 1);
for (std::size_t i = 0; i < nrows; i++) {
mat(i, 0) = phreeqc_mat.ids[i];
for (std::size_t j = 0; j < ncols; ++j) {
mat(i, j + 1) = phreeqc_mat.values[i][j];
}
}
Rcpp::colnames(mat) = Rcpp::wrap(phreeqc_mat.names);
return mat;
}
static inline Rcpp::List matToGrid(RInside &R, const Rcpp::NumericMatrix &mat,
const Rcpp::NumericMatrix &grid) {
Rcpp::Function pqc_to_grid_R("pqc_to_grid");
return pqc_to_grid_R(mat, grid);
}
static inline std::map<int, std::string>
replaceRawKeywordIDs(std::map<int, std::string> raws) {
for (auto &raw : raws) {
std::string &s = raw.second;
// find at beginning of line '*_RAW' followed by a number and change this
// number to 1
std::regex re(R"((RAW\s+)(\d+))");
s = std::regex_replace(s, re, "RAW 1");
}
return raws;
}
static inline uint32_t getSolutionCount(std::unique_ptr<PhreeqcInit> &phreeqc,
const Rcpp::List &initial_grid) {
PhreeqcInit::ModulesArray mod_array;
Rcpp::Function unique_R("unique");
std::vector<int> row_ids =
Rcpp::as<std::vector<int>>(unique_R(initial_grid["ID"]));
// std::vector<std::uint32_t> sizes_vec(sizes.begin(), sizes.end());
// Rcpp::Function modify_sizes("modify_module_sizes");
// sizes_vec = Rcpp::as<std::vector<std::uint32_t>>(
// modify_sizes(sizes_vec, phreeqc_mat, initial_grid));
// std::copy(sizes_vec.begin(), sizes_vec.end(), sizes.begin());
return phreeqc->getModuleSizes(row_ids)[POET_SOL];
}
static std::string readFile(const std::string &path) {
std::string string_rpath(PATH_MAX, '\0');
if (realpath(path.c_str(), string_rpath.data()) == nullptr) {
throw std::runtime_error("Failed to resolve the realpath to file " + path);
}
std::ifstream file(string_rpath);
if (!file.is_open()) {
throw std::runtime_error("Failed to open file: " + path);
}
std::stringstream buffer;
buffer << file.rdbuf();
return buffer.str();
}
void InitialList::prepareGrid(const Rcpp::List &grid_input) {
// parse input values
Rcpp::Function unique_R("unique");
std::string script;
std::string database;
if (grid_input.containsElementNamed(
GRID_MEMBER_STR(GridMembers::PQC_SCRIPT_FILE))) {
script = readFile(Rcpp::as<std::string>(
grid_input[GRID_MEMBER_STR(GridMembers::PQC_SCRIPT_FILE)]));
} else {
script = Rcpp::as<std::string>(
grid_input[GRID_MEMBER_STR(GridMembers::PQC_SCRIPT_STRING)]);
}
if (grid_input.containsElementNamed(
GRID_MEMBER_STR(GridMembers::PQC_DB_FILE))) {
database = readFile(Rcpp::as<std::string>(
grid_input[GRID_MEMBER_STR(GridMembers::PQC_DB_FILE)]));
} else {
database = Rcpp::as<std::string>(
grid_input[GRID_MEMBER_STR(GridMembers::PQC_DB_STRING)]);
}
this->database = database;
Rcpp::NumericMatrix grid_def =
grid_input[GRID_MEMBER_STR(GridMembers::GRID_DEF)];
Rcpp::NumericVector grid_size =
grid_input[GRID_MEMBER_STR(GridMembers::GRID_SIZE)];
// Rcpp::NumericVector constant_cells = grid["constant_cells"].;
this->n_rows = grid_def.nrow();
this->n_cols = grid_def.ncol();
this->s_cols = grid_size[0];
this->s_rows = grid_size[1];
this->dim = n_cols == 1 ? 1 : 2;
if (this->n_cols > 1 && this->n_rows == 1) {
throw std::runtime_error(
"Dimensions of grid definition does not match the expected format "
"(n_rows > 1, n_cols >= 1).");
}
if (this->dim != grid_size.size()) {
throw std::runtime_error(
"Dimensions of grid does not match the dimension of grid definition.");
}
if (this->s_rows <= 0 || this->s_cols <= 0) {
throw std::runtime_error("Grid size must be positive.");
}
this->phreeqc = std::make_unique<PhreeqcInit>(database, script);
this->phreeqc_mat = pqcScriptToGrid(phreeqc, R);
this->initial_grid = matToGrid(R, this->phreeqc_mat, grid_def);
const uint32_t solution_count = getSolutionCount(phreeqc, this->initial_grid);
std::vector<std::string> colnames =
Rcpp::as<std::vector<std::string>>(this->initial_grid.names());
this->transport_names = std::vector<std::string>(
colnames.begin() + 1,
colnames.begin() + 1 + solution_count); // skip ID
std::map<int, std::string> pqc_raw_dumps;
pqc_raw_dumps = replaceRawKeywordIDs(phreeqc->raw_dumps());
this->pqc_ids =
Rcpp::as<std::vector<int>>(unique_R(this->initial_grid["ID"]));
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));
}
}
} // namespace poet

135
src/Init/InitialList.cpp Normal file
View File

@ -0,0 +1,135 @@
#include "InitialList.hpp"
#include "DataStructures/NamedVector.hpp"
#include <Rcpp/internal/wrap.h>
#include <Rcpp/iostream/Rstreambuf.h>
#include <Rcpp/proxy/ProtectedProxy.h>
#include <Rcpp/vector/instantiation.h>
#include <cstdint>
#include <string>
#include <vector>
namespace poet {
void InitialList::initializeFromList(const Rcpp::List &setup) {
prepareGrid(setup[grid_key]);
initDiffusion(setup[diffusion_key]);
initChemistry(setup[chemistry_key]);
}
void InitialList::importList(const Rcpp::List &setup, bool minimal) {
this->dim = Rcpp::as<uint8_t>(setup[static_cast<int>(ExportList::GRID_DIM)]);
Rcpp::NumericVector grid_specs =
setup[static_cast<int>(ExportList::GRID_SPECS)];
this->n_rows = static_cast<std::uint32_t>(grid_specs[0]);
this->n_cols = static_cast<std::uint32_t>(grid_specs[1]);
Rcpp::NumericVector spatial =
setup[static_cast<int>(ExportList::GRID_SPATIAL)];
this->s_rows = spatial[0];
this->s_cols = spatial[1];
this->constant_cells = Rcpp::as<std::vector<uint32_t>>(
setup[static_cast<int>(ExportList::GRID_CONSTANT)]);
this->porosity = Rcpp::as<std::vector<double>>(
setup[static_cast<int>(ExportList::GRID_POROSITY)]);
if (!minimal) {
this->initial_grid =
Rcpp::List(setup[static_cast<int>(ExportList::GRID_INITIAL)]);
}
this->transport_names = Rcpp::as<std::vector<std::string>>(
setup[static_cast<int>(ExportList::DIFFU_TRANSPORT)]);
if (!minimal) {
this->boundaries =
Rcpp::List(setup[static_cast<int>(ExportList::DIFFU_BOUNDARIES)]);
this->inner_boundaries =
Rcpp::List(setup[static_cast<int>(ExportList::DIFFU_INNER_BOUNDARIES)]);
this->alpha_x =
Rcpp::List(setup[static_cast<int>(ExportList::DIFFU_ALPHA_X)]);
this->alpha_y =
Rcpp::List(setup[static_cast<int>(ExportList::DIFFU_ALPHA_Y)]);
}
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_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)]);
this->interp_species = Rcpp::as<NamedVector<uint32_t>>(
setup[static_cast<int>(ExportList::CHEM_INTERP_SPECIES)]);
this->chem_hooks =
Rcpp::as<Rcpp::List>(setup[static_cast<int>(ExportList::CHEM_HOOKS)]);
}
Rcpp::List InitialList::exportList() {
Rcpp::List out(static_cast<int>(ExportList::ENUM_SIZE));
out[static_cast<int>(ExportList::GRID_DIM)] = this->dim;
out[static_cast<int>(ExportList::GRID_SPECS)] =
Rcpp::NumericVector::create(this->n_rows, this->n_cols);
out[static_cast<int>(ExportList::GRID_SPATIAL)] =
Rcpp::NumericVector::create(this->s_rows, this->s_cols);
out[static_cast<int>(ExportList::GRID_CONSTANT)] =
Rcpp::wrap(this->constant_cells);
out[static_cast<int>(ExportList::GRID_POROSITY)] = Rcpp::wrap(this->porosity);
out[static_cast<int>(ExportList::GRID_INITIAL)] = this->initial_grid;
out[static_cast<int>(ExportList::DIFFU_TRANSPORT)] =
Rcpp::wrap(this->transport_names);
out[static_cast<int>(ExportList::DIFFU_BOUNDARIES)] = this->boundaries;
out[static_cast<int>(ExportList::DIFFU_INNER_BOUNDARIES)] =
this->inner_boundaries;
out[static_cast<int>(ExportList::DIFFU_ALPHA_X)] = this->alpha_x;
out[static_cast<int>(ExportList::DIFFU_ALPHA_Y)] = this->alpha_y;
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_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);
out[static_cast<int>(ExportList::CHEM_HOOKS)] = this->chem_hooks;
return out;
}
} // namespace poet

240
src/Init/InitialList.hpp Normal file
View File

@ -0,0 +1,240 @@
#pragma once
#include "Base/RInsidePOET.hpp"
#include "DataStructures/NamedVector.hpp"
#include "POETInit.hpp"
#include <Rcpp/vector/instantiation.h>
#include <set>
#include <tug/Boundary.hpp>
#include <RInside.h>
#include <Rcpp.h>
#include <array>
#include <cstddef>
#include <cstdint>
#include <PhreeqcInit.hpp>
#include <map>
#include <memory>
#include <string>
#include <utility>
#include <vector>
#include <DataStructures/Field.hpp>
namespace poet {
using TugType = double;
class InitialList {
public:
InitialList(RInside &R) : R(R){};
void initializeFromList(const Rcpp::List &setup);
void importList(const Rcpp::List &setup, bool minimal = false);
Rcpp::List exportList();
Field getInitialGrid() const { return Field(this->initial_grid); }
private:
RInside &R;
enum class ExportList {
GRID_DIM,
GRID_SPECS,
GRID_SPATIAL,
GRID_CONSTANT,
GRID_POROSITY,
GRID_INITIAL,
DIFFU_TRANSPORT,
DIFFU_BOUNDARIES,
DIFFU_INNER_BOUNDARIES,
DIFFU_ALPHA_X,
DIFFU_ALPHA_Y,
CHEM_DATABASE,
CHEM_FIELD_HEADER,
CHEM_PQC_SCRIPTS,
CHEM_PQC_IDS,
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,
ENUM_SIZE
};
// Grid members
static constexpr const char *grid_key = "Grid";
enum class GridMembers {
PQC_SCRIPT_STRING,
PQC_SCRIPT_FILE,
PQC_DB_STRING,
PQC_DB_FILE,
GRID_DEF,
GRID_SIZE,
CONSTANT_CELLS,
POROSITY,
ENUM_SIZE
};
static constexpr std::size_t size_GridMembers =
static_cast<std::size_t>(InitialList::GridMembers::ENUM_SIZE);
static constexpr std::array<const char *, size_GridMembers>
GridMembersString = {"pqc_in_string", "pqc_in_file", "pqc_db_string",
"pqc_db_file", "grid_def", "grid_size",
"constant_cells", "porosity"};
constexpr const char *GRID_MEMBER_STR(GridMembers member) const {
return GridMembersString[static_cast<std::size_t>(member)];
}
std::unique_ptr<PhreeqcInit> phreeqc;
void prepareGrid(const Rcpp::List &grid_input);
std::uint8_t dim{0};
std::uint32_t n_cols{0};
std::uint32_t n_rows{0};
double s_cols{0};
double s_rows{0};
std::vector<std::uint32_t> constant_cells;
std::vector<double> porosity;
Rcpp::List initial_grid;
// No export
Rcpp::NumericMatrix phreeqc_mat;
public:
struct DiffusionInit {
using BoundaryMap =
std::map<std::string, std::vector<tug::BoundaryElement<TugType>>>;
using InnerBoundaryMap =
std::map<std::string,
std::map<std::pair<std::uint32_t, std::uint32_t>, TugType>>;
uint8_t dim;
std::uint32_t n_cols;
std::uint32_t n_rows;
double s_cols;
double s_rows;
std::vector<std::uint32_t> constant_cells;
std::vector<std::string> transport_names;
BoundaryMap boundaries;
InnerBoundaryMap inner_boundaries;
Field alpha_x;
Field alpha_y;
};
DiffusionInit getDiffusionInit() const;
private:
// Diffusion members
static constexpr const char *diffusion_key = "Diffusion";
enum class DiffusionMembers {
BOUNDARIES,
INNER_BOUNDARIES,
ALPHA_X,
ALPHA_Y,
ENUM_SIZE
};
static constexpr std::size_t size_DiffusionMembers =
static_cast<std::size_t>(InitialList::DiffusionMembers::ENUM_SIZE);
static constexpr std::array<const char *, size_DiffusionMembers>
DiffusionMembersString = {"boundaries", "inner_boundaries", "alpha_x",
"alpha_y"};
constexpr const char *DIFFU_MEMBER_STR(DiffusionMembers member) const {
return DiffusionMembersString[static_cast<std::size_t>(member)];
}
void initDiffusion(const Rcpp::List &diffusion_input);
std::pair<Rcpp::List, Rcpp::List>
resolveBoundaries(const Rcpp::List &boundaries_list,
const Rcpp::List &inner_boundaries);
Rcpp::List boundaries;
Rcpp::List inner_boundaries;
Rcpp::List alpha_x;
Rcpp::List alpha_y;
std::vector<std::string> transport_names;
// Chemistry Members
static constexpr const char *chemistry_key = "Chemistry";
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_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;
NamedVector<std::uint32_t> interp_species;
Rcpp::List chem_hooks;
const std::set<std::string> hook_name_list{"dht_fill", "dht_fuzz",
"interp_pre", "interp_post"};
public:
struct ChemistryHookFunctions {
RHookFunction<bool> dht_fill;
RHookFunction<std::vector<double>> dht_fuzz;
RHookFunction<std::vector<std::size_t>> interp_pre;
RHookFunction<bool> interp_post;
};
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::map<int, POETConfig> pqc_config;
// std::vector<std::string> pqc_sol_order;
NamedVector<std::uint32_t> dht_species;
NamedVector<std::uint32_t> interp_species;
ChemistryHookFunctions hooks;
};
ChemistryInit getChemistryInit() const;
};
} // namespace poet

228
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,103 @@
*/
#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 "tug/Boundary.hpp"
#include "tug/Grid.hpp"
#include "tug/Simulation.hpp"
#include <array>
#include <cassert>
#include <mpi.h>
#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]);
if (!this->param_list.inner_boundaries[sol_name].empty()) {
boundary.setInnerBoundaries(this->param_list.inner_boundaries[sol_name]);
}
grid.setAlpha(alpha_x, alpha_y);
grid.setConcentrations(conc);
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; }

51
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,39 +69,23 @@ 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
*
*/
double transport_t = 0.f;
double transport_t = 0.;
};
} // namespace poet

83
src/initializer.cpp Normal file
View File

@ -0,0 +1,83 @@
#include "Init/InitialList.hpp"
#include "poet.hpp"
#include "Base/argh.hpp"
#include <cstdlib>
#include <RInside.h>
#include <Rcpp.h>
#include <iostream>
#include <string>
int main(int argc, char **argv) {
argh::parser cmdl({"-o", "--output"});
cmdl.parse(argc, argv);
if (cmdl[{"-h", "--help"}] || cmdl.pos_args().size() != 2) {
std::cout << "Usage: " << argv[0]
<< " [-o, --output output_file]"
" [-s, --setwd] "
" <script.R>"
<< std::endl;
return EXIT_SUCCESS;
}
RInside R(argc, argv);
R.parseEvalQ(init_r_library);
std::string input_script = cmdl.pos_args()[1];
std::string normalized_path_script;
std::string in_file_name;
std::string curr_path =
Rcpp::as<std::string>(Rcpp::Function("normalizePath")(Rcpp::wrap(".")));
try {
normalized_path_script =
Rcpp::as<std::string>(Rcpp::Function("normalizePath")(input_script));
in_file_name = Rcpp::as<std::string>(
Rcpp::Function("basename")(normalized_path_script));
std::cout << "Using script " << input_script << std::endl;
Rcpp::Function("source")(input_script);
} catch (Rcpp::exception &e) {
std::cerr << "Error while sourcing " << input_script << e.what()
<< std::endl;
return EXIT_FAILURE;
}
std::string output_file;
cmdl({"-o", "--output"},
curr_path + "/" +
in_file_name.substr(0, in_file_name.find_last_of('.')) + ".rds") >>
output_file;
if (cmdl[{"-s", "--setwd"}]) {
const std::string dir_path = Rcpp::as<std::string>(
Rcpp::Function("dirname")(normalized_path_script));
Rcpp::Function("setwd")(Rcpp::wrap(dir_path));
}
Rcpp::List setup = R["setup"];
poet::InitialList init(R);
init.initializeFromList(setup);
// replace file extension by .rds
Rcpp::Function save("saveRDS");
save(init.exportList(), Rcpp::wrap(output_file));
std::cout << "Saved result to " << output_file << std::endl;
// parseGrid(R, grid, results);
return EXIT_SUCCESS;
}

573
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,123 +20,237 @@
** 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 "DataStructures/Field.hpp"
#include "Init/InitialList.hpp"
#include "Transport/DiffusionModule.hpp"
#include <poet.hpp>
#include <RInside.h>
#include <Rcpp.h>
#include <cstdint>
#include <Rcpp/DataFrame.h>
#include <Rcpp/Function.h>
#include <Rcpp/vector/instantiation.h>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <string>
#include <vector>
#include <memory>
#include <mpi.h>
#include <optional>
#include <string>
#include "Base/argh.hpp"
#include <poet.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});
}
static std::unique_ptr<Rcpp::List> global_rt_setup;
return out_map;
// we need some layz evaluation, as we can't define the functions before the R
// runtime is initialized
static std::optional<Rcpp::Function> master_init_R;
static std::optional<Rcpp::Function> master_iteration_end_R;
static std::optional<Rcpp::Function> store_setup_R;
static void init_global_functions(RInside &R) {
R.parseEval(kin_r_library);
master_init_R = Rcpp::Function("master_init");
master_iteration_end_R = Rcpp::Function("master_iteration_end");
store_setup_R = Rcpp::Function("StoreSetup");
}
// 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");
// R["TMP"] = chem.asSEXP();
// R.parseEval("mysetup$state_C <- TMP");
// }
enum ParseRet { PARSER_OK, PARSER_ERROR, PARSER_HELP };
ParseRet parseInitValues(char **argv, RuntimeParameters &params) {
// initialize argh object
argh::parser cmdl(argv);
// if user asked for help
if (cmdl[{"help", "h"}]) {
if (MY_RANK == 0) {
MSG("Todo");
MSG("See README.md for further information.");
}
return ParseRet::PARSER_HELP;
}
// if positional arguments are missing
if (!cmdl(3)) {
if (MY_RANK == 0) {
ERRMSG("POET needs 3 positional arguments: ");
ERRMSG("1) the R script defining your simulation runtime.");
ERRMSG("2) the initial .rds file generated by poet_init.");
ERRMSG("3) the directory prefix where to save results and profiling");
}
return ParseRet::PARSER_ERROR;
}
// 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;
}
}
// 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;
}
}
params.print_progressbar = cmdl[{"P", "progress"}];
/*Parse work package size*/
cmdl("work-package-size", CHEM_DEFAULT_WORK_PACKAGE_SIZE) >>
params.work_package_size;
/* Parse DHT arguments */
params.use_dht = cmdl["dht"];
params.use_interp = cmdl["interp"];
// cout << "CPP: DHT is " << ( dht_enabled ? "ON" : "OFF" ) << '\n';
cmdl("dht-size", CHEM_DHT_SIZE_PER_PROCESS_MB) >> params.dht_size;
// cout << "CPP: DHT size per process (Byte) = " << dht_size_per_process <<
// endl;
cmdl("dht-snaps", 0) >> params.dht_snaps;
params.use_interp = cmdl["interp"];
cmdl("interp-size", 100) >> params.interp_size;
cmdl("interp-min", 5) >> params.interp_min_entries;
cmdl("interp-bucket-entries", 20) >> params.interp_bucket_entries;
if (MY_RANK == 0) {
// MSG("Complete results storage is " + BOOL_PRINT(simparams.store_result));
MSG("Work Package Size: " + std::to_string(params.work_package_size));
MSG("DHT is " + BOOL_PRINT(params.use_dht));
if (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(params.dht_size));
MSG("DHT save snapshots is " + BOOL_PRINT(params.dht_snaps));
// MSG("DHT load file is " + chem_params.dht_file);
}
if (params.use_interp) {
MSG("PHT interpolation enabled: " + BOOL_PRINT(params.use_interp));
MSG("PHT interp-size = " + std::to_string(params.interp_size));
MSG("PHT interp-min = " + std::to_string(params.interp_min_entries));
MSG("PHT interp-bucket-entries = " +
std::to_string(params.interp_bucket_entries));
}
}
std::string init_file;
std::string runtime_file;
cmdl(1) >> runtime_file;
cmdl(2) >> init_file;
cmdl(3) >> params.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_;
global_rt_setup = std::make_unique<Rcpp::List>();
*global_rt_setup = source(runtime_file, Rcpp::Named("local", true));
*global_rt_setup = global_rt_setup->operator[]("value");
params.timesteps =
Rcpp::as<std::vector<double>>(global_rt_setup->operator[]("timesteps"));
} catch (const std::exception &e) {
ERRMSG("Error while parsing R scripts: " + std::string(e.what()));
return ParseRet::PARSER_ERROR;
}
return ParseRet::PARSER_OK;
}
// 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(
"mysetup$state_T <- setNames(data.frame(matrix(TMP, nrow=" +
std::to_string(trans.GetRequestedVecSize()) + ")), TMP_PROPS)"));
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(
"mysetup$state_C <- setNames(data.frame(matrix(TMP, nrow=" +
std::to_string(chem.GetRequestedVecSize()) + ")), TMP_PROPS)"));
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"];
}
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);
static Rcpp::List RunMasterLoop(RInsidePOET &R, const RuntimeParameters &params,
DiffusionModule &diffusion,
ChemistryModule &chem) {
// 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);
}
inline double RunMasterLoop(SimParams &params, RInside &R,
const GridParams &g_params, uint32_t nxyz_master) {
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);
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());
if (params.getNumParams().print_progressbar) {
if (params.print_progressbar) {
chem.setProgressBarPrintout(true);
}
@ -143,47 +259,33 @@ inline double RunMasterLoop(SimParams &params, RInside &R,
double dSimTime{0};
for (uint32_t iter = 1; iter < maxiter + 1; iter++) {
double start_t = MPI_Wtime();
uint32_t tick = 0;
// 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) + "]"));
const 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");
/* 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();
writeFieldsToR(R, diffusion.getField(), chem.GetField());
diffusion.getField().update(chem.GetField());
R["req_dt"] = dt;
R["simtime"] = (sim_time += dt);
R.parseEval("mysetup$req_dt <- req_dt");
R.parseEval("mysetup$simtime <- simtime");
chem.simulate(dt);
// MDL master_iteration_end just writes on disk state_T and
// 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)");
call_master_iter_end(R, diffusion.getField(), chem.getField());
diffusion.getField().update(chem.getField());
MSG("End of *coupling* iteration " + std::to_string(iter) + "/" +
std::to_string(maxiter));
@ -194,175 +296,146 @@ inline double RunMasterLoop(SimParams &params, RInside &R,
dSimTime += end_t - start_t;
} // END SIMULATION LOOP
R.parseEvalQ("profiling <- list()");
Rcpp::List chem_profiling;
chem_profiling["simtime"] = chem.GetChemistryTime();
chem_profiling["loop"] = chem.GetMasterLoopTime();
chem_profiling["sequential"] = chem.GetMasterSequentialTime();
chem_profiling["idle_master"] = chem.GetMasterIdleTime();
chem_profiling["idle_worker"] = Rcpp::wrap(chem.GetWorkerIdleTimings());
chem_profiling["phreeqc_time"] = Rcpp::wrap(chem.GetWorkerPhreeqcTimings());
R["simtime_chemistry"] = chem.GetChemistryTime();
R.parseEvalQ("profiling$simtime_chemistry <- simtime_chemistry");
Rcpp::List diffusion_profiling;
diffusion_profiling["simtime"] = diffusion.getTransportTime();
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["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["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["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.use_dht) {
chem_profiling["dht_hits"] = Rcpp::wrap(chem.GetWorkerDHTHits());
chem_profiling["dht_evictions"] = Rcpp::wrap(chem.GetWorkerDHTEvictions());
chem_profiling["dht_get_time"] = Rcpp::wrap(chem.GetWorkerDHTGetTimings());
chem_profiling["dht_fill_time"] =
Rcpp::wrap(chem.GetWorkerDHTFillTimings());
}
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"] =
if (params.use_interp) {
chem_profiling["interp_w"] =
Rcpp::wrap(chem.GetWorkerInterpolationWriteTimings());
chem_profiling["interp_r"] =
Rcpp::wrap(chem.GetWorkerInterpolationReadTimings());
chem_profiling["interp_g"] =
Rcpp::wrap(chem.GetWorkerInterpolationGatherTimings());
chem_profiling["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_profiling["interp_calls"] =
Rcpp::wrap(chem.GetWorkerInterpolationCalls());
chem_profiling["interp_cached"] = Rcpp::wrap(chem.GetWorkerPHTCacheHits());
}
Rcpp::List profiling;
profiling["simtime"] = dSimTime;
profiling["chemistry"] = chem_profiling;
profiling["diffusion"] = diffusion_profiling;
chem.MasterLoopBreak();
diffusion.end();
return dSimTime;
return profiling;
}
int main(int argc, char *argv[]) {
double dSimTime, sim_end;
int world_size, world_rank;
int world_size;
MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD, &world_size);
{
MPI_Comm_size(MPI_COMM_WORLD, &world_size);
MPI_Comm_rank(MPI_COMM_WORLD, &MY_RANK);
MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);
RInsidePOET &R = RInsidePOET::getInstance();
RInsidePOET &R = RInsidePOET::getInstance();
if (world_rank == 0) {
MSG("Running POET version " + std::string(poet_version));
}
if (world_rank > 0) {
{
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;
}
// 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();
if (MY_RANK == 0) {
MSG("Running POET version " + std::string(poet_version));
}
RuntimeParameters run_params;
switch (parseInitValues(argv, run_params)) {
case ParseRet::PARSER_ERROR:
case ParseRet::PARSER_HELP:
MPI_Finalize();
return 0;
case ParseRet::PARSER_OK:
break;
}
InitialList init_list(R);
init_list.importList(run_params.init_params, MY_RANK != 0);
MSG("RInside initialized on process " + std::to_string(MY_RANK));
std::cout << std::flush;
MPI_Barrier(MPI_COMM_WORLD);
MSG("finished, cleanup of process " + std::to_string(world_rank));
ChemistryModule chemistry(run_params.work_package_size,
init_list.getChemistryInit(), MPI_COMM_WORLD);
MPI_Finalize();
const ChemistryModule::SurrogateSetup surr_setup = {
init_list.getInitialGrid().GetProps(),
run_params.use_dht,
run_params.dht_size,
run_params.use_interp,
run_params.interp_bucket_entries,
run_params.interp_size,
run_params.interp_min_entries};
return EXIT_SUCCESS;
chemistry.masterEnableSurrogates(surr_setup);
if (MY_RANK > 0) {
chemistry.WorkerLoop();
} else {
init_global_functions(R);
// R.parseEvalQ("mysetup <- setup");
// // if (MY_RANK == 0) { // get timestep vector from
// // grid_init function ... //
*global_rt_setup =
master_init_R.value()(*global_rt_setup, run_params.out_dir,
init_list.getInitialGrid().asSEXP());
// MDL: store all parameters
// MSG("Calling R Function to store calling parameters");
// R.parseEvalQ("StoreSetup(setup=mysetup)");
MSG("Init done on process with rank " + std::to_string(MY_RANK));
// MPI_Barrier(MPI_COMM_WORLD);
DiffusionModule diffusion(init_list.getDiffusionInit(),
init_list.getInitialGrid());
chemistry.masterSetField(init_list.getInitialGrid());
Rcpp::List profiling = RunMasterLoop(R, run_params, diffusion, chemistry);
MSG("finished simulation loop");
R["profiling"] = profiling;
R["setup"] = *global_rt_setup;
string r_vis_code;
r_vis_code =
"saveRDS(profiling, file=paste0(setup$out_dir,'/timings.rds'));";
R.parseEval(r_vis_code);
MSG("Done! Results are stored as R objects into <" + run_params.out_dir +
"/timings.rds>");
}
}
/*Loading Dependencies*/
// TODO: kann raus
R.parseEvalQ(kin_r_library);
MSG("finished, cleanup of process " + std::to_string(MY_RANK));
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;
}
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.parseEval(master_init_code);
GridParams g_params(R);
params.initVectorParams(R);
// MDL: store all parameters
if (world_rank == 0) {
MSG("Calling R Function to store calling parameters");
R.parseEvalQ("StoreSetup(setup=mysetup)");
}
if (world_rank == 0) {
MSG("Init done on process with rank " + std::to_string(world_rank));
}
// MPI_Barrier(MPI_COMM_WORLD);
uint32_t nxyz_master = (world_size == 1 ? g_params.total_n : 1);
dSimTime = RunMasterLoop(params, R, g_params, nxyz_master);
MSG("finished simulation loop");
MSG("start timing profiling");
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);
MSG("Done! Results are stored as R objects into <" + params.getOutDir() +
"/timings.rds>");
MPI_Barrier(MPI_COMM_WORLD);
MSG("finished, cleanup of process " + std::to_string(world_rank));
MPI_Finalize();
if (world_rank == 0) {
if (MY_RANK == 0) {
MSG("done, bye!");
}

93
src/poet.hpp.in
View File

@ -1,10 +1,95 @@
#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 <cstdint>
#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@)";
#endif // POET_H
static const inline std::string init_r_library = R"(@R_INIT_LIB@)";
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"};
constexpr uint32_t CHEM_DEFAULT_WORK_PACKAGE_SIZE = 32;
constexpr uint32_t CHEM_DHT_SIZE_PER_PROCESS_MB = 1.5E3;
struct RuntimeParameters {
std::string out_dir;
std::vector<double> timesteps;
bool print_progressbar;
uint32_t work_package_size;
Rcpp::List init_params;
bool use_dht;
std::uint32_t dht_size;
std::uint8_t dht_snaps;
bool use_interp;
std::uint32_t interp_size;
std::uint32_t interp_min_entries;
std::uint32_t interp_bucket_entries;
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);
};
};

4
test/CMakeLists.txt
View File

@ -2,8 +2,10 @@ file(GLOB test_SRC
CONFIGURE_DEPENDS
"*.cpp" "*.c")
find_package(doctest REQUIRED)
add_executable(testPOET ${test_SRC})
target_link_libraries(testPOET doctest poetlib)
target_link_libraries(testPOET doctest::doctest POETLib)
target_include_directories(testPOET PRIVATE "${PROJECT_SOURCE_DIR}/src")
get_filename_component(TEST_RInsideSourceFile "RInsidePOET_funcs.R" REALPATH)

3
test/RInsidePOET_funcs.R
View File

@ -13,10 +13,13 @@ bool_named_vec <- function(input) {
}
simple_field <- function(field) {
field <- as.data.frame(field, check.names = FALSE)
field$Na <- 0
return(field)
}
extended_field <- function(field, additional) {
field <- as.data.frame(field, check.names = FALSE)
additional <- as.data.frame(additional, check.names = FALSE)
return(field + additional)
}

5
test/testField.cpp
View File

@ -1,15 +1,12 @@
#include <Rcpp.h>
#include <algorithm>
#include <array>
#include <cstddef>
#include <cstdint>
#include <doctest/doctest.h>
#include <iostream>
#include <string>
#include <vector>
#include <Base/RInsidePOET.hpp>
#include <DataStructures/DataStructures.hpp>
#include <DataStructures/Field.hpp>
#include "testDataStructures.hpp"

3
test/testNamedVector.cpp
View File

@ -1,11 +1,10 @@
#include <Rcpp.h>
#include <cstddef>
#include <doctest/doctest.h>
#include <utility>
#include <vector>
#include <Base/RInsidePOET.hpp>
#include <DataStructures/DataStructures.hpp>
#include <DataStructures/NamedVector.hpp>
#include "testDataStructures.hpp"

11
util/ci.Dockerfile
View File

@ -1,11 +0,0 @@
FROM debian
RUN apt-get update && export DEBIAN_FRONTEND=noninteractive \
&& apt-get install -y \
cmake \
git \
libeigen3-dev \
libopenmpi-dev \
r-cran-rcpp \
r-cran-rinside \
libssl-dev

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