max/poet
1
0
Fork 0
mirror of https://git.gfz-potsdam.de/naaice/poet.git synced 2025-12-16 12:54:50 +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 74bafca4b8
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 \ RUN sudo apt-get update && export DEBIAN_FRONTEND=noninteractive \
&& sudo apt-get install -y \ && sudo apt-get install -y \
cmake-curses-gui \ cmake \
clangd \
git \ git \
libeigen3-dev \ libeigen3-dev \
libopenmpi-dev \ libopenmpi-dev \
r-cran-rcpp \ r-base-dev
r-cran-rinside
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/ build/
/.cache/ /.cache/
.vscode .vscode
.codechecker

9
.gitmodules vendored
View File

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

12
CMakeLists.txt
View File

@ -1,12 +1,11 @@
# Version 3.9+ offers new MPI package variables cmake_minimum_required(VERSION 3.14)
cmake_minimum_required(VERSION 3.9)
project(POET project(POET
LANGUAGES CXX C LANGUAGES CXX C
DESCRIPTION "A coupled reactive transport simulator") DESCRIPTION "A coupled reactive transport simulator")
# specify the C++ standard # specify the C++ standard
set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD 20)
set(CMAKE_CXX_STANDARD_REQUIRED True) set(CMAKE_CXX_STANDARD_REQUIRED True)
set(CMAKE_INSTALL_RPATH_USE_LINK_PATH 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() get_poet_version()
# set(GCC_CXX_FLAGS "-D STRICT_R_HEADERS") add_definitions(${GCC_CXX_FLAGS})
find_package(MPI REQUIRED) find_package(MPI REQUIRED)
find_package(RRuntime REQUIRED) find_package(RRuntime REQUIRED)
@ -29,15 +26,14 @@ add_subdirectory(bench)
set(TUG_ENABLE_TESTING OFF CACHE BOOL "" FORCE) set(TUG_ENABLE_TESTING OFF CACHE BOOL "" FORCE)
add_subdirectory(ext/tug EXCLUDE_FROM_ALL) 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) option(POET_ENABLE_TESTING "Build test suite for POET" OFF)
if (POET_ENABLE_TESTING) if (POET_ENABLE_TESTING)
add_subdirectory(ext/doctest EXCLUDE_FROM_ALL)
add_subdirectory(test) add_subdirectory(test)
endif() endif()
option(BUILD_DOC "Build documentation with doxygen" OFF) 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: The following external header library is shipped with POET:
- **argh** - https://github.com/adishavit/argh (BSD license) - **argh** - https://github.com/adishavit/argh (BSD license)
- **PhreeqcRM** with patches from GFZ - - **IPhreeqc** with patches from GFZ -
https://www.usgs.gov/software/phreeqc-version-3 - https://github.com/usgs-coupled/iphreeqc -
https://git.gfz-potsdam.de/mluebke/phreeqcrm-gfz https://git.gfz-potsdam.de/naaice/iphreeqc
- **tug** - https://git.gfz-potsdam.de/sec34/tug - **tug** - https://git.gfz-potsdam.de/naaice/tug
## Installation ## Installation
@ -35,6 +35,7 @@ To compile POET you need several software to be installed:
- MPI-Implementation (tested with OpenMPI and MVAPICH) - MPI-Implementation (tested with OpenMPI and MVAPICH)
- R language and environment - R language and environment
- CMake 3.9+ - CMake 3.9+
- Eigen3 3.4+ (required by `tug`)
- *optional*: `doxygen` with `dot` bindings for documentiation - *optional*: `doxygen` with `dot` bindings for documentiation
The following R libraries must then be installed, which will get the 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 ```sh
poet poet
├── bin ├── bin
│ └── poet │   ├── poet
├── R_lib │   └── poet_init
│ └── kin_r_library.R
└── share └── share
└── poet └── poet
└── bench ├── barite
├── barite │   ├── barite_200.rds
│ ├── barite_interp_eval.R │   ├── barite_200_rt.R
│ ├── barite.pqi │   ├── barite_het.rds
│ ├── barite.R │   └── barite_het_rt.R
│ └── db_barite.dat ├── dolo
├── dolo │   ├── dolo_inner_large.rds
│ ├── dolo_diffu_inner_large.R │   ├── dolo_inner_large_rt.R
│ ├── dolo_diffu_inner.R │   ├── dolo_interp.rds
│ ├── dolo_inner.pqi │   └── dolo_interp_rt.R
│ ├── dolo_interp_long.R └── surfex
│ └── phreeqc_kin.dat ├── PoetEGU_surfex_500.rds
└── surfex └── PoetEGU_surfex_500_rt.R
├── ExBase.pqi
├── ex.R
├── SMILE_2021_11_01_TH.dat
├── SurfExBase.pqi
└── surfex.R
``` ```
The R libraries will be loaded at runtime and the paths are hardcoded With the installation of POET, two executables are provided:
absolute paths inside `poet.cpp`. So, if you consider to move - `poet` - the main executable to run simulations
`bin/poet` either change paths of the R source files and recompile - `poet_init` - a preprocessor to generate input files for POET from R scripts
POET or also move `R_lib/*` relative to the binary.
The benchmarks consist of input scripts, which are provided as .R files. Preprocessed benchmarks can be found in the `share/poet` directory with an
Additionally, Phreeqc scripts and their corresponding databases are required, according *runtime* setup. More on those files and how to create them later.
stored as .pqi and .dat files, respectively.
## Running ## Running
Run POET by `mpirun ./poet <OPTIONS> <SIMFILE> <OUTPUT_DIRECTORY>` Run POET by `mpirun ./poet [OPTIONS] <RUNFILE> <SIMFILE> <OUTPUT_DIRECTORY>`
where: where:
- **OPTIONS** - runtime parameters (explained below) - **OPTIONS** - POET options (explained below)
- **SIMFILE** - simulation described as R script (e.g. - **RUNFILE** - Runtime parameters described as R script
`<POET_INSTALL_DIR>/share/poet/bench/dolo/dolo_interp_long.R`) - **SIMFILE** - Simulation input prepared by `poet_init`
- **OUTPUT_DIRECTORY** - path, where all output of POET should be stored - **OUTPUT_DIRECTORY** - path, where all output of POET should be stored
### Runtime options ### POET options
The following parameters can be set: The following parameters can be set:
| Option | Value | Description | | Option | Value | Description |
|-----------------------------|--------------|--------------------------------------------------------------------------------------------------------------------------| |-----------------------------|--------------|--------------------------------------------------------------------------------------------------------------------------|
| **--work-package-size=** | _1..n_ | size of work packages (defaults to _5_) | | **--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** | | enabling DHT usage (defaults to _OFF_) |
| **--dht-strategy=** | _0-1_ | change DHT strategy. **NOT IMPLEMENTED YET** (Defaults to _0_) | | **--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_) | | **--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 ### Example: Running from scratch
We will continue the above example and start a simulation with We will continue the above example and start a simulation with *barite_het*,
`dolo_diffu_inner.R`. As transport a simple fixed-coefficient diffusion is used. which simulation files can be found in
It's a 2D, 100x100 grid, simulating 10 time steps. To start the simulation with `<INSTALL_DIR>/share/poet/barite/barite_het*`. As transport a heterogeneous
4 processes `cd` into your previously installed POET-dir diffusion is used. It's a small 2D grid, 2x5 grid, simulating 50 time steps with
`<POET_INSTALL_DIR>/bin` and run: 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 ```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 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: resulting call would look like this:
```sh ```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() ## About the usage of MPI_Wtime()
Implemented time measurement functions uses `MPI_Wtime()`. Some 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 ### this program; if not, write to the Free Software Foundation, Inc., 51
### Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. ### Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
master_init <- function(setup, out_dir, init_field) {
## 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")
## Setup the directory where we will store the results ## Setup the directory where we will store the results
verb <- FALSE if (!dir.exists(out_dir)) {
if (local_rank == 0) { dir.create(out_dir)
verb <- TRUE ## verbosity loading MUFITS results msgm("created directory ", out_dir)
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)
}
} else { } 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$iter <- 1
setup$maxiter <- setup$iterations
setup$timesteps <- setup$timesteps setup$timesteps <- setup$timesteps
setup$maxiter <- length(setup$timesteps)
setup$iterations <- setup$maxiter
setup$simulation_time <- 0 setup$simulation_time <- 0
if (is.null(setup[["store_result"]])) { if (is.null(setup[["store_result"]])) {
setup$store_result <- TRUE setup$store_result <- TRUE
} }
if (setup$store_result) { 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"]])) { if (is.null(setup[["out_save"]])) {
setup$out_save <- seq(1, setup$iterations) setup$out_save <- seq(1, setup$iterations)
} }
} }
setup$out_dir <- out_dir
return(setup) return(setup)
} }
## This function, called only by master, stores on disk the last ## This function, called only by master, stores on disk the last
## calculated time step if store_result is TRUE and increments the ## calculated time step if store_result is TRUE and increments the
## iteration counter ## iteration counter
master_iteration_end <- function(setup) { master_iteration_end <- function(setup, state_T, state_C) {
iter <- setup$iter iter <- setup$iter
# print(iter)
## max digits for iterations ## max digits for iterations
dgts <- as.integer(ceiling(log10(setup$iterations + 1))) dgts <- as.integer(ceiling(log10(setup$maxiter)))
## string format to use in sprintf ## string format to use in sprintf
fmt <- paste0("%0", dgts, "d") fmt <- paste0("%0", dgts, "d")
## Write on disk state_T and state_C after every iteration ## Write on disk state_T and state_C after every iteration
## comprised in setup$out_save ## comprised in setup$out_save
if (setup$store_result) { if (setup$store_result) {
if (iter %in% setup$out_save) { if (iter %in% setup$out_save) {
nameout <- paste0(fileout, "/iter_", sprintf(fmt=fmt, iter), ".rds") nameout <- paste0(setup$out_dir, "/iter_", sprintf(fmt = fmt, iter), ".rds")
info <- list(
tr_req_dt = as.integer(setup$req_dt) state_T <- data.frame(state_T, check.names = FALSE)
## tr_allow_dt = setup$allowed_dt, state_C <- data.frame(state_C, check.names = FALSE)
## tr_inniter = as.integer(setup$inniter)
)
saveRDS(list( saveRDS(list(
T = setup$state_T, C = setup$state_C, T = state_T,
simtime = as.integer(setup$simtime), C = state_C,
tr_info = info simtime = as.integer(setup$simulation_time)
), file = nameout) ), file = nameout)
msgm("results stored in <", 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 setup$iter <- setup$iter + 1
return(setup) 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 ## Attach the name of the calling function to the message displayed on
## R's stdout ## R's stdout
msgm <- function(...) { msgm <- function(...) {
if (local_rank == 0) { prefix <- paste0("R: ")
fname <- as.list(sys.call(-1))[[1]] cat(paste(prefix, ..., "\n"))
prefix <- paste0("R: ", fname, " ::")
cat(paste(prefix, ..., "\n"))
}
invisible() invisible()
} }
## Function called by master R process to store on disk all relevant ## Function called by master R process to store on disk all relevant
## parameters for the simulation ## parameters for the simulation
StoreSetup <- function(setup) { StoreSetup <- function(setup, filesim, out_dir) {
to_store <- vector(mode = "list", length = 4) to_store <- vector(mode = "list", length = 4)
## names(to_store) <- c("Sim", "Flow", "Transport", "Chemistry", "DHT") ## names(to_store) <- c("Sim", "Flow", "Transport", "Chemistry", "DHT")
names(to_store) <- c("Sim", "Transport", "DHT", "Cmdline") names(to_store) <- c("Sim", "Transport", "DHT", "Cmdline")
## read the setup R file, which is sourced in kin.cpp ## read the setup R file, which is sourced in kin.cpp
tmpbuff <- file(filesim, "r") tmpbuff <- file(filesim, "r")
setupfile <- readLines(tmpbuff) setupfile <- readLines(tmpbuff)
close.connection(tmpbuff) close.connection(tmpbuff)
to_store$Sim <- setupfile to_store$Sim <- setupfile
## to_store$Flow <- list( ## to_store$Flow <- list(
## snapshots = setup$snapshots, ## snapshots = setup$snapshots,
## gridfile = setup$gridfile, ## gridfile = setup$gridfile,
@ -282,22 +148,22 @@ StoreSetup <- function(setup) {
to_store$DHT <- FALSE to_store$DHT <- FALSE
} }
if (dht_enabled) { if (dht_enabled) {
to_store$DHT <- list( to_store$DHT <- list(
enabled = dht_enabled, enabled = dht_enabled,
log = dht_log log = dht_log
#signif = dht_final_signif, # signif = dht_final_signif,
#proptype = dht_final_proptype # proptype = dht_final_proptype
) )
} else { } else {
to_store$DHT <- FALSE to_store$DHT <- FALSE
} }
saveRDS(to_store, file = paste0(fileout, "/setup.rds")) saveRDS(to_store, file = paste0(fileout, "/setup.rds"))
msgm("initialization stored in ", paste0(fileout, "/setup.rds")) msgm("initialization stored in ", paste0(fileout, "/setup.rds"))
} }
GetWorkPackageSizesVector <- function(n_packages, package_size, len) { 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))) 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(dolo)
add_subdirectory(surfex) add_subdirectory(surfex)
add_subdirectory(barite)

26
bench/barite/CMakeLists.txt
View File

@ -1,8 +1,20 @@
install(FILES # Create a list of files
barite.R set(bench_files
barite_interp_eval.R barite_200.R
barite.pqi barite_het.R
db_barite.dat
DESTINATION
share/poet/bench/barite
) )
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 * List of Files
- =barite.R=: POET input script for a 20x20 simulation grid - =barite_het.R=: POET input script with homogeneous zones for a 5x2 simulation
- =barite_interp_eval.R=: POET input script for a 400x200 simulation grid
- =barite_200.R=: POET input script for a 200x200 simulation
grid grid
- =db_barite.dat=: PHREEQC database containing the kinetic expressions - =db_barite.dat=: PHREEQC database containing the kinetic expressions
for barite and celestite, stripped down from =phreeqc.dat= 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 SOLUTION 1
units mol/kgw units mol/kgw
water 1 water 1
temperature 25 temperature 25
pH 7 pH 7
pe 10.799 PURE 1
Ba 0.1 Celestite 0.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
END 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") s_cols <- 1
input_script <- normalizePath("../share/poet/bench/barite/barite.pqi") s_rows <- 1
## database <- normalizePath("/home/work/simR/Rphree/poetsims/Sims/Hans/db_barite.dat")
## input_script <- normalizePath("/home/work/simR/Rphree/poetsims/Sims/Hans/barite.pqi")
################################################################# grid_def <- matrix(2, nrow = rows, ncol = cols)
## Section 1 ##
## Grid initialization ##
#################################################################
n <- 200 # Define grid configuration for POET model
m <- 200 grid_setup <- list(
pqc_in_file = "./barite.pqi",
init_cell <- list( pqc_db_file = "./db_barite.dat", # Path to the database file for Phreeqc
"H" = 110.0124, grid_def = grid_def, # Definition of the grid, containing IDs according to the Phreeqc input script
"O" = 55.5087, grid_size = c(s_rows, s_cols), # Size of the grid in meters
"Charge" = -1.216307845207e-09, constant_cells = c() # IDs of cells with constant concentration
"Ba" = 1.E-12,
"Cl" = 2.E-12,
"S(6)" = 6.204727095976e-04,
"Sr" = 6.204727095976e-04,
"Barite" = 0.001,
"Celestite" = 1
) )
grid <- list( bound_length <- 2
n_cells = c(n, m),
s_cells = c(1, 1), bound_def <- list(
type = "scratch", "type" = rep("constant", bound_length),
init_cell = as.data.frame(init_cell, check.names = FALSE), "sol_id" = rep(3, bound_length),
props = names(init_cell), "cell" = seq(1, bound_length)
database = database,
input_script = input_script
) )
homogenous_alpha <- 1e-6
################################################################## diffusion_setup <- list(
## Section 2 ## boundaries = list(
## Diffusion parameters and boundary conditions ## "W" = bound_def,
################################################################## "N" = bound_def
),
## initial conditions alpha_x = homogenous_alpha,
alpha_y = homogenous_alpha
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) ##
#################################################################
## DHT significant digits
dht_species <- c( dht_species <- c(
"H" = 7, "H" = 7,
"O" = 7, "O" = 7,
@ -112,27 +46,13 @@ dht_species <- c(
"Celestite" = 4 "Celestite" = 4
) )
chemistry <- list( chemistry_setup <- list(
database = database, dht_species = dht_species
input_script = input_script,
dht_species = dht_species
) )
################################################################# # Define a setup list for simulation configuration
## Section 4 ##
## Putting all those things together ##
#################################################################
iterations <- 50
dt <- 100
setup <- list( setup <- list(
grid = grid, Grid = grid_setup, # Parameters related to the grid structure
diffusion = diffusion, Diffusion = diffusion_setup, # Parameters related to the diffusion process
chemistry = chemistry, Chemistry = chemistry_setup
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))
) )

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 set(bench_files
dolo_diffu_inner.R dolo_inner_large.R
dolo_diffu_inner_large.R dolo_interp.R
dolo_inner.pqi
dolo_interp_long.R
phreeqc_kin.dat
DESTINATION
share/poet/bench/dolo
) )
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 * List of Files
- =dolo_diffu_inner.R=: POET input script for a 100x100 simulation - =dolo_interp.R=: POET input script for a 400x200 simulation
grid
- =dolo_interp_long.R=: POET input script for a 400x200 simulation
grid grid
- =dolo_diffu_inner_large.R=: POET input script for a 400x200
simulation grid
- =phreeqc_kin.dat=: PHREEQC database containing the kinetic expressions - =phreeqc_kin.dat=: PHREEQC database containing the kinetic expressions
for dolomite and celestite, stripped down from =phreeqc.dat= for dolomite and celestite, stripped down from =phreeqc.dat=
- =dol.pqi=: PHREEQC input script for the chemical system - =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 * 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 SOLUTION 1
units mol/kgw units mol/kgw
temp 25.0 water 1
water 1 temperature 25
pH 9.91 charge pH 7
pe 4.0 pe 4
C 1.2279E-04
Ca 1.2279E-04
Cl 1E-12
Mg 1E-12
PURE 1 PURE 1
O2g -0.1675 10 Calcite 0.0 1
KINETICS 1 END
Calcite
-m 0.000207 RUN_CELLS
-parms 0.0032 -cells 1
-tol 1e-10
Dolomite COPY solution 1 2
-m 0.0
-parms 0.00032 PURE 2
-tol 1e-10 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 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 set(bench_files
ExBase.pqi # surfex.R
ex.R # ex.R
surfex.R PoetEGU_surfex_500.R
SurfExBase.pqi
SMILE_2021_11_01_TH.dat
DESTINATION
share/poet/bench/surfex
) )
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 Z 0.0012585
Y 0.0009418 Y 0.0009418
END 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 - =ex.R=: POET input script for a 100x100 simulation grid, only
exchange exchange
- =ExBase.pqi=: PHREEQC input script for the =ex.R= model - =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 considering both cation exchange and surface complexation
- =SurfExBase.pqi=: PHREEQC input script for the =surfex.R= model - =SurfExBase.pqi=: PHREEQC input script for the =surfex.R= model
- =SMILE_2021_11_01_TH.dat=: PHREEQC database containing the - =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 Z 0.0012585
Y 0.0009418 Y 0.0009418
END 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") grid_def <- matrix(1, nrow = rows, ncol = cols)
input_script <- normalizePath("./ExBase.pqi")
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
)
################################################################# bound_def <- list(
## Section 1 ## "type" = rep("constant", cols),
## Grid initialization ## "sol_id" = rep(2, cols),
################################################################# "cell" = seq(1, cols)
)
n <- 100 diffusion_setup <- list(
m <- 100 boundaries = list(
"N" = bound_def
types <- c("scratch", "phreeqc", "rds") ),
alpha_x = 1e-6,
init_cell <- list(H = 1.476571028625e-01, alpha_y = 1e-6
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"
) )
################################################################## chemistry_setup <- list()
## 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
# Define a setup list for simulation configuration
setup <- list( setup <- list(
grid = grid, Grid = grid_setup, # Parameters related to the grid structure
diffusion = diffusion, Diffusion = diffusion_setup, # Parameters related to the diffusion process
chemistry = chemistry, Chemistry = chemistry_setup # Parameters related to the chemistry process
iterations = iterations, )
timesteps = rep(dt, iterations),
store_result = TRUE
)

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") grid_def <- matrix(1, nrow = rows, ncol = cols)
input_script <- normalizePath("../share/poet/bench/surfex/SurfExBase.pqi")
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
)
################################################################# bound_def <- list(
## Section 1 ## "type" = rep("constant", cols),
## Grid initialization ## "sol_id" = rep(2, cols),
################################################################# "cell" = seq(1, cols)
)
n <- 1000 diffusion_setup <- list(
m <- 1000 boundaries = list(
"N" = bound_def
types <- c("scratch", "phreeqc", "rds") ),
alpha_x = 1e-6,
init_cell <- list(H = 1.476571028625e-01, alpha_y = 1e-6
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"
) )
################################################################## chemistry_setup <- list()
## 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
# Define a setup list for simulation configuration
setup <- list( setup <- list(
grid = grid, Grid = grid_setup, # Parameters related to the grid structure
diffusion = diffusion, Diffusion = diffusion_setup, # Parameters related to the diffusion process
chemistry = chemistry, Chemistry = chemistry_setup # Parameters related to the chemistry process
iterations = iterations, )
timesteps = rep(dt, iterations),
store_result = TRUE,
out_save = c(5, iterations)
)

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 doxygen_add_docs(doxygen
${PROJECT_SOURCE_DIR}/src ${PROJECT_SOURCE_DIR}/src
${PROJECT_SOURCE_DIR}/README.md ${PROJECT_SOURCE_DIR}/README.md
${PROJECT_SOURCE_DIR}/docs/Input_Scripts.md
${PROJECT_SOURCE_DIR}/docs/Output.md ${PROJECT_SOURCE_DIR}/docs/Output.md
COMMENT "Generate html pages") COMMENT "Generate html pages")
endif() 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 and possible to read out within a R runtime with
`readRDS("timings.rds")`. There you will find the following values: `readRDS("timings.rds")`. There you will find the following values:
| Value | Description | | Value | Description |
|--------------------|----------------------------------------------------------------------------| | --------- | -------------------------------------------------------------------------- |
| simtime | time spent in whole simulation loop without any initialization and cleanup | | simtime | time spent in whole simulation loop without any initialization and cleanup |
| simtime\_transport | measured time in *transport* subroutine | | chemistry | measured time in *chemistry* subroutine |
| simtime\_chemistry | measured time in *chemistry* subroutine (actual parallelized part) | | diffusion | measured time in *diffusion* subroutine |
### chemistry subsetting ### Chemistry subsetting
If running parallel there are also measured timings which are subsets of | Value | Description |
*simtime\_​chemistry*. | ------------- | --------------------------------------------------------- |
| 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 | #### DHT usage
|-----------------------|-----------------------------------------------------------|
| 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}
If running in parallel and with activated DHT, two more timings and also If running in parallel and with activated DHT, two more timings and also
some profiling about the DHT usage are given: some profiling about the DHT usage are given:
| Value | Description | | Value | Description |
|-----------------|---------------------------------------------------------| | --------------- | ------------------------------------------------------- |
| dht\_fill\_time | time to write data to DHT | | dht\_hits | count of data points retrieved from 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\_evictions | count of data points evicted by another write operation | | 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 <RInside.h>
#include <Rcpp.h> #include <Rcpp.h>
#include <Rinternals.h>
#include <cstddef> #include <cstddef>
#include <exception> #include <exception>
#include <optional> #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 { template <typename... Args> T operator()(Args... args) const {
if (func.has_value()) { if (func.has_value()) {
return (Rcpp::as<T>(this->func.value()(args...))); 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(); } bool isValid() const { return this->func.has_value(); }
SEXP asSEXP() const { return Rcpp::as<SEXP>(this->func.value()); }
private: private:
std::optional<Rcpp::Function> func; 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 add_library(POETLib)
Base/Grid.cpp
Base/SimParams.cpp target_sources(POETLib
Chemistry/ChemistryModule.cpp PRIVATE
Chemistry/MasterFunctions.cpp Init/InitialList.cpp
Chemistry/WorkerFunctions.cpp Init/GridInit.cpp
Chemistry/SurrogateModels/DHT_Wrapper.cpp Init/DiffusionInit.cpp
Chemistry/SurrogateModels/DHT.c Init/ChemistryInit.cpp
Chemistry/SurrogateModels/HashFunctions.cpp
Chemistry/SurrogateModels/InterpolationModule.cpp
Chemistry/SurrogateModels/ProximityHashTable.cpp
DataStructures/Field.cpp DataStructures/Field.cpp
Transport/DiffusionModule.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 target_include_directories(POETLib PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}")
MPI::MPI_C target_link_libraries(
${MATH_LIBRARY} POETLib
RRuntime PUBLIC RRuntime
PhreeqcRM PUBLIC IPhreeqcPOET
tug 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) # target_link_libraries(poetlib PUBLIC
set(PHREEQCRM_DISABLE_OPENMP ON CACHE BOOL "" FORCE) # 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) # mark_as_advanced(PHREEQCRM_BUILD_MPI PHREEQCRM_DISABLE_OPENMP)
target_compile_definitions(poetlib PRIVATE DHT_STATISTICS) # set(PHREEQCRM_DISABLE_OPENMP ON CACHE BOOL "" FORCE)
endif()
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) # if(POET_DHT_DEBUG)
target_compile_definitions(poetlib PRIVATE POET_PHT_ADD) # target_compile_definitions(poetlib PRIVATE DHT_STATISTICS)
endif() # 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/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) configure_file(poet.hpp.in poet.hpp @ONLY)
add_executable(poet poet.cpp) 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}") 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 "ChemistryModule.hpp"
#include "PhreeqcEngine.hpp"
#include "SurrogateModels/DHT_Wrapper.hpp" #include "SurrogateModels/DHT_Wrapper.hpp"
#include "SurrogateModels/Interpolation.hpp" #include "SurrogateModels/Interpolation.hpp"
#include <PhreeqcRM.h>
#include <algorithm> #include <algorithm>
#include <cassert> #include <cassert>
#include <cstddef>
#include <cstdint> #include <cstdint>
#include <map>
#include <memory> #include <memory>
#include <mpi.h> #include <mpi.h>
#include <stdexcept>
#include <string> #include <string>
#include <utility>
#include <vector> #include <vector>
constexpr uint32_t MB_FACTOR = 1E6;
std::vector<double> std::vector<double>
inverseDistanceWeighting(const std::vector<std::int32_t> &to_calc, inverseDistanceWeighting(const std::vector<std::int32_t> &to_calc,
const std::vector<double> &from, 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) { const std::vector<std::vector<double>> &output) {
std::vector<double> results = from; std::vector<double> results = from;
const std::uint32_t buffer_size = input.size() + 1; // const std::uint32_t buffer_size = input.size() + 1;
double buffer[buffer_size]; // double buffer[buffer_size];
double from_rescaled; // double from_rescaled;
const std::uint32_t data_set_n = input.size(); const std::uint32_t data_set_n = input.size();
double rescaled[to_calc.size()][data_set_n + 1]; double rescaled[to_calc.size()][data_set_n + 1];
double weights[data_set_n]; double weights[data_set_n];
// rescaling over all key elements // 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]; const auto output_comp_i = to_calc[key_comp_i];
// rescale input between 0 and 1 // 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]; 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], const double max = *std::max_element(rescaled[key_comp_i],
rescaled[key_comp_i] + data_set_n + 1); 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] = rescaled[key_comp_i][point_i] =
((max - min) != 0 ((max - min) != 0
? (rescaled[key_comp_i][point_i] - min) / (max - min) ? (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 // calculate distances for each data set
double inv_sum = 0; 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; 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( distance += std::pow(
rescaled[key_comp_i][point_i] - rescaled[key_comp_i][data_set_n], 2); 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_h = false;
// bool has_o = 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]; const auto output_comp_i = to_calc[key_comp_i];
double key_delta = 0; double key_delta = 0;
@ -91,7 +87,7 @@ inverseDistanceWeighting(const std::vector<std::int32_t> &to_calc,
// has_o = true; // 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]; key_delta += weights[j] * output[j][output_comp_i];
} }
@ -158,25 +154,30 @@ inverseDistanceWeighting(const std::vector<std::int32_t> &to_calc,
return results; return results;
} }
poet::ChemistryModule::ChemistryModule(uint32_t nxyz, uint32_t wp_size, poet::ChemistryModule::ChemistryModule(
std::uint32_t maxiter, uint32_t wp_size_, const InitialList::ChemistryInit chem_params,
const ChemistryParams &chem_param, MPI_Comm communicator)
MPI_Comm communicator) : group_comm(communicator), wp_size(wp_size_), params(chem_params) {
: PhreeqcRM(nxyz, 1), group_comm(communicator), wp_size(wp_size), MPI_Comm_rank(communicator, &comm_rank);
params(chem_param) { MPI_Comm_size(communicator, &comm_size);
MPI_Comm_size(communicator, &this->comm_size); this->is_sequential = comm_size == 1;
MPI_Comm_rank(communicator, &this->comm_rank); this->is_master = comm_rank == 0;
this->is_sequential = (this->comm_size == 1); this->n_cells = chem_params.total_grid_cells;
this->is_master = (this->comm_rank == 0);
if (!is_sequential && is_master) { if (!is_master) {
MPI_Bcast(&wp_size, 1, MPI_UINT32_T, 0, this->group_comm); for (const auto &[pqc_id, pqc_config] : chem_params.pqc_config) {
MPI_Bcast(&maxiter, 1, MPI_UINT32_T, 0, this->group_comm); 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() { 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( void poet::ChemistryModule::initializeDHT(
uint32_t size_mb, const NamedVector<std::uint32_t> &key_species) { uint32_t size_mb, const NamedVector<std::uint32_t> &key_species) {
constexpr uint32_t MB_FACTOR = 1E6; constexpr uint32_t MB_FACTOR = 1E6;
@ -368,7 +207,7 @@ void poet::ChemistryModule::initializeDHT(
if (key_species.empty()) { if (key_species.empty()) {
std::vector<std::uint32_t> default_signif( 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); 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->dht = new DHT_Wrapper(dht_comm, dht_size, map_copy, key_indices,
this->prop_names, params.hooks, 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)); 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; std::vector<int32_t> species_indices;
species_indices.reserve(key_species.size()); species_indices.reserve(key_species.size());
const auto test = key_species.getNames();
for (const auto &species : key_species.getNames()) { for (const auto &species : key_species.getNames()) {
auto it = std::find(to_compare.begin(), to_compare.end(), species); auto it = std::find(to_compare.begin(), to_compare.end(), species);
if (it == to_compare.end()) { if (it == to_compare.end()) {
@ -469,9 +310,9 @@ void poet::ChemistryModule::initializeInterp(
if (key_species.empty()) { if (key_species.empty()) {
map_copy = this->dht->getKeySpecies(); 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 = 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 ? InterpolationModule::COARSE_DIFF
: 0); : 0);
map_copy[i] = signif; map_copy[i] = signif;

120
src/Chemistry/ChemistryModule.hpp
View File

@ -2,16 +2,17 @@
#ifndef CHEMISTRYMODULE_H_ #ifndef CHEMISTRYMODULE_H_
#define CHEMISTRYMODULE_H_ #define CHEMISTRYMODULE_H_
#include "../Base/SimParams.hpp" #include "DataStructures/Field.hpp"
#include "../DataStructures/DataStructures.hpp" #include "DataStructures/NamedVector.hpp"
#include "ChemistryDefs.hpp"
#include "Init/InitialList.hpp"
#include "SurrogateModels/DHT_Wrapper.hpp" #include "SurrogateModels/DHT_Wrapper.hpp"
#include "SurrogateModels/Interpolation.hpp" #include "SurrogateModels/Interpolation.hpp"
#include <IrmResult.h> #include "PhreeqcEngine.hpp"
#include <PhreeqcRM.h>
#include <array> #include <array>
#include <cstddef>
#include <cstdint> #include <cstdint>
#include <map> #include <map>
#include <memory> #include <memory>
@ -24,7 +25,7 @@ namespace poet {
* \brief Wrapper around PhreeqcRM to provide POET specific parallelization with * \brief Wrapper around PhreeqcRM to provide POET specific parallelization with
* easy access. * easy access.
*/ */
class ChemistryModule : public PhreeqcRM { class ChemistryModule {
public: public:
/** /**
* Creates a new instance of Chemistry module with given grid cell count, work * 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 wp_size Count of grid cells to fill each work package at maximum.
* \param communicator MPI communicator to distribute work in. * \param communicator MPI communicator to distribute work in.
*/ */
ChemistryModule(uint32_t nxyz, uint32_t wp_size, std::uint32_t maxiter, ChemistryModule(uint32_t wp_size,
const ChemistryParams &chem_param, MPI_Comm communicator); const InitialList::ChemistryInit chem_params,
MPI_Comm communicator);
/** /**
* Deconstructor, which frees DHT data structure if used. * Deconstructor, which frees DHT data structure if used.
*/ */
~ChemistryModule(); ~ChemistryModule();
/** void masterSetField(Field field);
* 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);
/** /**
* Run the chemical simulation with parameters set. * 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 * Returns all known species names, including not only aqueous species, but
* also equilibrium, exchange, surface and kinetic reactants. * 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. * Return the accumulated runtime in seconds for chemical simulation.
*/ */
auto GetChemistryTime() const { return this->chem_t; } auto GetChemistryTime() const { return this->chem_t; }
/** void setFilePadding(std::uint32_t maxiter) {
* Create a new worker instance inside given MPI communicator. this->file_pad =
* static_cast<std::uint8_t>(std::ceil(std::log10(maxiter + 1)));
* 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);
/** struct SurrogateSetup {
* Default work package size. std::vector<std::string> prop_names;
*/
static constexpr uint32_t CHEM_DEFAULT_WP_SIZE = 5; 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 * Intended to alias input parameters for grid initialization with a single
@ -118,17 +125,6 @@ public:
DHT_SNAPS_ITEREND //!< output snapshots after each iteration 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. * **Only called by workers!** Start the worker listening loop.
*/ */
@ -243,9 +239,7 @@ protected:
const NamedVector<std::uint32_t> &key_species); const NamedVector<std::uint32_t> &key_species);
enum { enum {
CHEM_INIT, CHEM_FIELD_INIT,
CHEM_WP_SIZE,
CHEM_INIT_SPECIES,
CHEM_DHT_ENABLE, CHEM_DHT_ENABLE,
CHEM_DHT_SIGNIF_VEC, CHEM_DHT_SIGNIF_VEC,
CHEM_DHT_SNAPS, CHEM_DHT_SNAPS,
@ -294,7 +288,7 @@ protected:
using worker_list_t = std::vector<struct worker_info_s>; using worker_list_t = std::vector<struct worker_info_s>;
using workpointer_t = std::vector<double>::iterator; using workpointer_t = std::vector<double>::iterator;
void MasterRunParallel(); void MasterRunParallel(double dt);
void MasterRunSequential(); void MasterRunSequential();
void MasterSendPkgs(worker_list_t &w_list, workpointer_t &work_pointer, 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 WorkerPerfToMaster(int type, const struct worker_s &timings);
void WorkerMetricsToMaster(int type); void WorkerMetricsToMaster(int type);
IRM_RESULT WorkerRunWorkPackage(WorkPackage &work_package, double dSimTime, void WorkerRunWorkPackage(WorkPackage &work_package, double dSimTime,
double dTimestep); double dTimestep);
std::vector<uint32_t> CalculateWPSizesVector(uint32_t n_cells, std::vector<uint32_t> CalculateWPSizesVector(uint32_t n_cells,
uint32_t wp_size) const; uint32_t wp_size) const;
@ -344,7 +338,7 @@ protected:
bool is_sequential; bool is_sequential;
bool is_master; bool is_master;
uint32_t wp_size{CHEM_DEFAULT_WP_SIZE}; uint32_t wp_size;
bool dht_enabled{false}; bool dht_enabled{false};
int dht_snaps_type{DHT_SNAPS_DISABLED}; int dht_snaps_type{DHT_SNAPS_DISABLED};
std::string dht_file_out_dir; std::string dht_file_out_dir;
@ -372,7 +366,7 @@ protected:
bool print_progessbar{false}; bool print_progessbar{false};
std::uint32_t file_pad; std::uint8_t file_pad{1};
double chem_t{0.}; double chem_t{0.};
@ -382,11 +376,9 @@ protected:
Field chem_field; Field chem_field;
static constexpr int MODULE_COUNT = 5; const InitialList::ChemistryInit params;
const ChemistryParams &params; std::map<int, std::unique_ptr<PhreeqcEngine>> phreeqc_instances;
std::array<std::uint32_t, MODULE_COUNT> speciesPerModule{};
}; };
} // namespace poet } // namespace poet

59
src/Chemistry/MasterFunctions.cpp
View File

@ -1,12 +1,9 @@
#include "ChemistryModule.hpp" #include "ChemistryModule.hpp"
#include <PhreeqcRM.h>
#include <algorithm> #include <algorithm>
#include <cstddef> #include <cstddef>
#include <cstdint> #include <cstdint>
#include <mpi.h> #include <mpi.h>
#include <stdexcept>
#include <vector> #include <vector>
std::vector<uint32_t> 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()}; double start_t{MPI_Wtime()};
if (this->is_sequential) { if (this->is_sequential) {
MasterRunSequential(); MasterRunSequential();
return; return;
} }
MasterRunParallel(); MasterRunParallel(dt);
double end_t{MPI_Wtime()}; double end_t{MPI_Wtime()};
this->chem_t += end_t - start_t; this->chem_t += end_t - start_t;
} }
void poet::ChemistryModule::MasterRunSequential() { void poet::ChemistryModule::MasterRunSequential() {
std::vector<double> shuffled_field = // std::vector<double> shuffled_field =
shuffleField(chem_field.AsVector(), n_cells, prop_count, 1); // shuffleField(chem_field.AsVector(), n_cells, prop_count, 1);
std::vector<std::vector<double>> input; // std::vector<std::vector<double>> input;
for (std::size_t i = 0; i < n_cells; i++) { // for (std::size_t i = 0; i < n_cells; i++) {
input.push_back( // input.push_back(
std::vector<double>(shuffled_field.begin() + (i * prop_count), // std::vector<double>(shuffled_field.begin() + (i * prop_count),
shuffled_field.begin() + ((i + 1) * prop_count))); // shuffled_field.begin() + ((i + 1) *
} // prop_count)));
// }
this->setDumpedField(input); // this->setDumpedField(input);
PhreeqcRM::RunCells(); // PhreeqcRM::RunCells();
this->getDumpedField(input); // this->getDumpedField(input);
shuffled_field.clear(); // shuffled_field.clear();
for (std::size_t i = 0; i < n_cells; i++) { // for (std::size_t i = 0; i < n_cells; i++) {
shuffled_field.insert(shuffled_field.end(), input[i].begin(), // shuffled_field.insert(shuffled_field.end(), input[i].begin(),
input[i].end()); // input[i].end());
} // }
std::vector<double> out_vec{shuffled_field}; // std::vector<double> out_vec{shuffled_field};
unshuffleField(shuffled_field, n_cells, prop_count, 1, out_vec); // unshuffleField(shuffled_field, n_cells, prop_count, 1, out_vec);
chem_field = out_vec; // chem_field = out_vec;
} }
void poet::ChemistryModule::MasterRunParallel() { void poet::ChemistryModule::MasterRunParallel(double dt) {
/* declare most of the needed variables here */ /* declare most of the needed variables here */
double seq_a, seq_b, seq_c, seq_d; double seq_a, seq_b, seq_c, seq_d;
double worker_chemistry_a, worker_chemistry_b; double worker_chemistry_a, worker_chemistry_b;
@ -360,7 +358,6 @@ void poet::ChemistryModule::MasterRunParallel() {
MPI_Barrier(this->group_comm); MPI_Barrier(this->group_comm);
double dt = this->PhreeqcRM::GetTimeStep();
static uint32_t iteration = 0; static uint32_t iteration = 0;
/* start time measurement of sequential part */ /* start time measurement of sequential part */
@ -466,3 +463,13 @@ poet::ChemistryModule::CalculateWPSizesVector(uint32_t n_cells,
return wp_sizes_vector; 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 "DHT_Wrapper.hpp"
#include "Init/InitialList.hpp"
#include "Rounding.hpp"
#include <algorithm> #include <algorithm>
#include <cassert> #include <cassert>
#include <cmath> #include <cmath>
#include <cstddef> #include <cstddef>
#include <cstdint> #include <cstdint>
#include <cstring> #include <cstring>
#include <iostream>
#include <stdexcept> #include <stdexcept>
#include <vector> #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 NamedVector<std::uint32_t> &key_species,
const std::vector<std::int32_t> &key_indices, const std::vector<std::int32_t> &key_indices,
const std::vector<std::string> &_output_names, const std::vector<std::string> &_output_names,
const ChemistryParams::Chem_Hook_Functions &_hooks, const InitialList::ChemistryHookFunctions &_hooks,
uint32_t data_count, bool _with_interp) uint32_t data_count, bool _with_interp)
: key_count(key_indices.size()), data_count(data_count), : key_count(key_indices.size()), data_count(data_count),
input_key_elements(key_indices), communicator(dht_comm), 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 #ifndef DHT_WRAPPER_H
#define DHT_WRAPPER_H #define DHT_WRAPPER_H
#include "../../Base/RInsidePOET.hpp" #include "Base/RInsidePOET.hpp"
#include "../../Base/SimParams.hpp" #include "DataStructures/NamedVector.hpp"
#include "../../DataStructures/DataStructures.hpp"
#include "../enums.hpp" #include "Chemistry/ChemistryDefs.hpp"
#include "HashFunctions.hpp"
#include "Init/InitialList.hpp"
#include "LookupKey.hpp" #include "LookupKey.hpp"
#include "Rounding.hpp"
#include <array> #include <array>
#include <cstdint> #include <cstdint>
#include <limits> #include <limits>
#include <string> #include <string>
#include <unordered_map>
#include <utility> #include <utility>
#include <vector> #include <vector>
@ -88,7 +87,7 @@ public:
const NamedVector<std::uint32_t> &key_species, const NamedVector<std::uint32_t> &key_species,
const std::vector<std::int32_t> &key_indices, const std::vector<std::int32_t> &key_indices,
const std::vector<std::string> &output_names, const std::vector<std::string> &output_names,
const ChemistryParams::Chem_Hook_Functions &hooks, const InitialList::ChemistryHookFunctions &hooks,
uint32_t data_count, bool with_interp); uint32_t data_count, bool with_interp);
/** /**
* @brief Destroy the dht wrapper object * @brief Destroy the dht wrapper object
@ -260,7 +259,7 @@ private:
const std::vector<std::string> &output_names; const std::vector<std::string> &output_names;
const ChemistryParams::Chem_Hook_Functions &hooks; const InitialList::ChemistryHookFunctions &hooks;
const bool with_interp; const bool with_interp;
DHT_ResultObject dht_results; DHT_ResultObject dht_results;

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

@ -3,14 +3,13 @@
#ifndef INTERPOLATION_H_ #ifndef INTERPOLATION_H_
#define INTERPOLATION_H_ #define INTERPOLATION_H_
#include "../../Base/SimParams.hpp" #include "DataStructures/NamedVector.hpp"
#include "../../DataStructures/DataStructures.hpp"
#include "DHT_Wrapper.hpp" #include "DHT_Wrapper.hpp"
#include "Init/InitialList.hpp"
#include "LookupKey.hpp" #include "LookupKey.hpp"
#include "Rounding.hpp" #include "Rounding.hpp"
#include <cassert>
#include <iostream>
#include <list> #include <list>
#include <memory> #include <memory>
#include <mpi.h> #include <mpi.h>
@ -22,7 +21,6 @@ extern "C" {
} }
#include <cstdint> #include <cstdint>
#include <functional>
#include <unordered_map> #include <unordered_map>
#include <vector> #include <vector>
@ -165,7 +163,7 @@ public:
const NamedVector<std::uint32_t> &interp_key_signifs, const NamedVector<std::uint32_t> &interp_key_signifs,
const std::vector<std::int32_t> &dht_key_indices, const std::vector<std::int32_t> &dht_key_indices,
const std::vector<std::string> &out_names, 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 }; enum result_status { RES_OK, INSUFFICIENT_DATA, NOT_NEEDED };
@ -261,7 +259,7 @@ private:
return out_key; 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> &out_names;
const std::vector<std::string> dht_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" // Time-stamp: "Last modified 2023-08-16 17:02:31 mluebke"
#include "Init/InitialList.hpp"
#include "Interpolation.hpp" #include "Interpolation.hpp"
#include "../../DataStructures/DataStructures.hpp"
#include "DHT_Wrapper.hpp" #include "DHT_Wrapper.hpp"
#include "DataStructures/NamedVector.hpp"
#include "HashFunctions.hpp" #include "HashFunctions.hpp"
#include "LookupKey.hpp" #include "LookupKey.hpp"
#include "Rounding.hpp" #include "Rounding.hpp"
@ -35,7 +36,7 @@ InterpolationModule::InterpolationModule(
const NamedVector<std::uint32_t> &interp_key_signifs, const NamedVector<std::uint32_t> &interp_key_signifs,
const std::vector<std::int32_t> &dht_key_indices, const std::vector<std::int32_t> &dht_key_indices,
const std::vector<std::string> &_out_names, 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), : dht_instance(dht), key_signifs(interp_key_signifs),
key_indices(dht_key_indices), min_entries_needed(min_entries_needed), key_indices(dht_key_indices), min_entries_needed(min_entries_needed),
dht_names(dht.getKeySpecies().getNames()), out_names(_out_names), 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/DHT_Wrapper.hpp"
#include "SurrogateModels/Interpolation.hpp" #include "SurrogateModels/Interpolation.hpp"
#include <IrmResult.h> #include "Chemistry/ChemistryDefs.hpp"
#include <algorithm> #include <algorithm>
#include <cassert>
#include <cmath>
#include <cstddef> #include <cstddef>
#include <cstdint> #include <cstdint>
#include <iomanip> #include <iomanip>
@ -17,6 +16,7 @@
#include <vector> #include <vector>
namespace poet { namespace poet {
inline std::string get_string(int root, MPI_Comm communicator) { inline std::string get_string(int root, MPI_Comm communicator) {
int count; int count;
MPI_Bcast(&count, 1, MPI_INT, root, communicator); MPI_Bcast(&count, 1, MPI_INT, root, communicator);
@ -45,13 +45,8 @@ void poet::ChemistryModule::WorkerLoop() {
PropagateFunctionType(func_type); PropagateFunctionType(func_type);
switch (func_type) { switch (func_type) {
case CHEM_INIT: { case CHEM_FIELD_INIT: {
RunInitFile(get_string(0, this->group_comm)); ChemBCast(&this->prop_count, 1, MPI_UINT32_T);
break;
}
case CHEM_INIT_SPECIES: {
Field dummy;
initializeField(dummy);
break; break;
} }
case CHEM_WORK_LOOP: { case CHEM_WORK_LOOP: {
@ -177,9 +172,7 @@ void poet::ChemistryModule::WorkerDoWork(MPI_Status &probe_status,
phreeqc_time_start = MPI_Wtime(); phreeqc_time_start = MPI_Wtime();
if (WorkerRunWorkPackage(s_curr_wp, current_sim_time, dt) != IRM_OK) { WorkerRunWorkPackage(s_curr_wp, current_sim_time, dt);
std::cerr << "Phreeqc error" << std::endl;
};
phreeqc_time_end = MPI_Wtime(); phreeqc_time_end = MPI_Wtime();
@ -285,42 +278,34 @@ void poet::ChemistryModule::WorkerReadDHTDump(
} }
} }
IRM_RESULT void poet::ChemistryModule::WorkerRunWorkPackage(WorkPackage &work_package,
poet::ChemistryModule::WorkerRunWorkPackage(WorkPackage &work_package, double dSimTime,
double dSimTime, double dTimestep) { double dTimestep) {
// check if we actually need to start phreeqc
std::vector<std::uint32_t> pqc_mapping;
for (std::size_t i = 0; i < work_package.size; i++) { for (std::size_t wp_id = 0; wp_id < work_package.size; wp_id++) {
if (work_package.mapping[i] == CHEM_PQC) { if (work_package.mapping[wp_id] != CHEM_PQC) {
pqc_mapping.push_back(i); 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, 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.h>
#include <Rcpp/DataFrame.h>
#include <Rcpp/exceptions.h>
#include <Rcpp/vector/instantiation.h>
#include <Rinternals.h> #include <Rinternals.h>
#include <cstddef> #include <cstddef>
#include <cstdint> #include <cstdint>
@ -117,40 +120,14 @@ poet::FieldColumn &poet::Field::operator[](const std::string &key) {
} }
SEXP poet::Field::asSEXP() const { SEXP poet::Field::asSEXP() const {
const std::size_t cols = this->props.size(); Rcpp::List output;
SEXP s_names = PROTECT(Rf_allocVector(STRSXP, cols)); for (const auto &elem : this->props) {
SEXP s_output = PROTECT(Rf_allocVector(VECSXP, cols)); const auto map_it = this->find(elem);
output[elem] = Rcpp::wrap(map_it->second);
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);
} }
SEXP s_rownames = PROTECT(Rf_allocVector(INTSXP, this->req_vec_size)); return output;
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;
} }
poet::Field &poet::Field::operator=(const FieldColumn &cont_field) { 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) { void poet::Field::fromSEXP(const SEXP &s_rhs) {
this->clear(); this->clear();
SEXP s_vec = PROTECT(Rf_coerceVector(s_rhs, VECSXP)); Rcpp::List in_list;
SEXP s_names = PROTECT(Rf_getAttrib(s_vec, R_NamesSymbol));
std::size_t cols = static_cast<std::size_t>(Rf_length(s_vec)); try {
in_list = Rcpp::List(s_rhs);
this->props.clear(); } catch (const Rcpp::exception &e) {
this->props.reserve(cols); throw std::runtime_error("Input cannot be casted as list.");
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});
} }
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_ #ifndef DATASTRUCTURES_H_
#define DATASTRUCTURES_H_ #define DATASTRUCTURES_H_
#include "../Chemistry/enums.hpp"
#include <Rcpp.h> #include <Rcpp.h>
#include <cassert>
#include <cinttypes>
#include <cstddef> #include <cstddef>
#include <iostream> #include <cstdint>
#include <list>
#include <string> #include <string>
#include <utility>
#include <vector> #include <vector>
namespace poet { 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>; 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) 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 ** 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 ** 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 ** Foundation; either version 2 of the License, or (at your option) any later
@ -19,174 +22,103 @@
*/ */
#include "DiffusionModule.hpp" #include "DiffusionModule.hpp"
#include "../Base/Macros.hpp"
#include <tug/BoundaryCondition.hpp> #include "Base/Macros.hpp"
#include <tug/Diffusion.hpp> #include "Init/InitialList.hpp"
#include <Rcpp.h> #include <Eigen/src/Core/Map.h>
#include <algorithm> #include <Eigen/src/Core/Matrix.h>
#include <cstdint> #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 <string>
#include <vector> #include <vector>
using namespace poet; 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 = { for (std::uint32_t i = 0; i < mat.cols(); i++) {
tug::bc::BC_SIDE_LEFT, tug::bc::BC_SIDE_RIGHT, tug::bc::BC_SIDE_TOP, for (std::uint32_t j = 0; j < mat.rows(); j++) {
tug::bc::BC_SIDE_BOTTOM}; vec[j * mat.cols() + i] = mat(j, i);
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);
} }
} }
this->t_field = Field(n_grid_cells, initial_values, prop_names); return vec;
// 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);
}
}
}
} }
void DiffusionModule::simulate(double dt) { static inline Eigen::MatrixX<TugType>
double sim_a_transport, sim_b_transport; 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(); for (std::uint32_t i = 0; i < n_cols; i++) {
for (std::uint32_t j = 0; j < n_rows; j++) {
MSG_NOENDL("DiffusionModule::simulate(): Starting diffusion ..."); mat(j, i) = vec[j * n_cols + i];
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());
} }
} }
t_field = field_2d; return mat;
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;
} }
void DiffusionModule::end() { // static constexpr double ZERO_MULTIPLICATOR = 10E-14;
// R["simtime_transport"] = transport_t;
// R.parseEvalQ("profiling$simtime_transport <- simtime_transport"); 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; } double DiffusionModule::getTransportTime() { return this->transport_t; }

51
src/Transport/DiffusionModule.hpp
View File

@ -21,20 +21,13 @@
#ifndef DIFFUSION_MODULE_H #ifndef DIFFUSION_MODULE_H
#define DIFFUSION_MODULE_H #define DIFFUSION_MODULE_H
#include "../Base/Grid.hpp" #include "DataStructures/Field.hpp"
#include "../Base/SimParams.hpp" #include "Init/InitialList.hpp"
#include "../DataStructures/DataStructures.hpp"
#include <tug/BoundaryCondition.hpp> #include <sys/types.h>
#include <tug/Diffusion.hpp>
#include <array>
#include <cmath>
#include <cstdint>
#include <string>
#include <vector>
namespace poet { namespace poet {
/** /**
* @brief Class describing transport simulation * @brief Class describing transport simulation
* *
@ -53,8 +46,8 @@ public:
* *
* @param R RRuntime object * @param R RRuntime object
*/ */
DiffusionModule(const poet::DiffusionParams &diffu_args, DiffusionModule(const InitialList::DiffusionInit &init_list, Field field)
const poet::GridParams &grid_params); : param_list(init_list), transport_field(field){};
/** /**
* @brief Run simulation for one iteration * @brief Run simulation for one iteration
@ -64,14 +57,6 @@ public:
*/ */
void simulate(double dt); void simulate(double dt);
/**
* @brief End simulation
*
* All measured timings are distributed to the R runtime
*
*/
void end();
/** /**
* @brief Get the transport time * @brief Get the transport time
* *
@ -84,39 +69,23 @@ public:
* *
* \return Reference to the diffusion field. * \return Reference to the diffusion field.
*/ */
Field &getField() { return this->t_field; } Field &getField() { return this->transport_field; }
private: private:
/** /**
* @brief Instance of RRuntime * @brief Instance of RRuntime
* *
*/ */
// RRuntime &R;
enum { DIM_1D = 1, DIM_2D }; InitialList::DiffusionInit param_list;
void initialize(const poet::DiffusionParams &args, Field transport_field;
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;
/** /**
* @brief time spent for transport * @brief time spent for transport
* *
*/ */
double transport_t = 0.f; double transport_t = 0.;
}; };
} // namespace poet } // 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) 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 ** 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 ** 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 ** 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. ** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/ */
#include "Base/Grid.hpp"
#include "Base/Macros.hpp" #include "Base/Macros.hpp"
#include "Base/RInsidePOET.hpp" #include "Base/RInsidePOET.hpp"
#include "Base/SimParams.hpp"
#include "Chemistry/ChemistryModule.hpp" #include "Chemistry/ChemistryModule.hpp"
#include "DataStructures/Field.hpp"
#include "Init/InitialList.hpp"
#include "Transport/DiffusionModule.hpp" #include "Transport/DiffusionModule.hpp"
#include <poet.hpp> #include <RInside.h>
#include <Rcpp.h> #include <Rcpp.h>
#include <cstdint> #include <Rcpp/DataFrame.h>
#include <Rcpp/Function.h>
#include <Rcpp/vector/instantiation.h>
#include <cstdlib> #include <cstdlib>
#include <memory>
#include <cstring>
#include <iostream>
#include <string>
#include <vector>
#include <mpi.h> #include <mpi.h>
#include <optional>
#include <string>
#include "Base/argh.hpp"
#include <poet.hpp>
using namespace std; using namespace std;
using namespace poet; using namespace poet;
using namespace Rcpp; using namespace Rcpp;
poet::ChemistryModule::SingleCMap DFToHashMap(const Rcpp::DataFrame &df) { static int MY_RANK = 0;
std::unordered_map<std::string, double> out_map;
vector<string> col_names = Rcpp::as<vector<string>>(df.names());
for (const auto &name : col_names) { static std::unique_ptr<Rcpp::List> global_rt_setup;
double val = df[name.c_str()];
out_map.insert({name, val});
}
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 // HACK: this is a step back as the order and also the count of fields is
// predefined, but it will change in the future // 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"] = Rcpp::wrap(trans.AsVector());
R["TMP_PROPS"] = Rcpp::wrap(trans.GetProps()); R["TMP_PROPS"] = Rcpp::wrap(trans.GetProps());
R.parseEval(std::string( R.parseEval(std::string("state_T <- setNames(data.frame(matrix(TMP, nrow=" +
"mysetup$state_T <- setNames(data.frame(matrix(TMP, nrow=" + std::to_string(trans.GetRequestedVecSize()) +
std::to_string(trans.GetRequestedVecSize()) + ")), TMP_PROPS)")); ")), TMP_PROPS)"));
R["TMP"] = Rcpp::wrap(chem.AsVector()); R["TMP"] = Rcpp::wrap(chem.AsVector());
R["TMP_PROPS"] = Rcpp::wrap(chem.GetProps()); R["TMP_PROPS"] = Rcpp::wrap(chem.GetProps());
R.parseEval(std::string( R.parseEval(std::string("state_C <- setNames(data.frame(matrix(TMP, nrow=" +
"mysetup$state_C <- setNames(data.frame(matrix(TMP, nrow=" + std::to_string(chem.GetRequestedVecSize()) +
std::to_string(chem.GetRequestedVecSize()) + ")), TMP_PROPS)")); ")), 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, static Rcpp::List RunMasterLoop(RInsidePOET &R, const RuntimeParameters &params,
const std::string &database_path) { DiffusionModule &diffusion,
chem.SetErrorHandlerMode(1); ChemistryModule &chem) {
chem.SetComponentH2O(false);
chem.SetRebalanceFraction(0.5);
chem.SetRebalanceByCell(true);
chem.UseSolutionDensityVolume(false);
chem.SetPartitionUZSolids(false);
// Set concentration units
// 1, mg/L; 2, mol/L; 3, kg/kgs
chem.SetUnitsSolution(2);
// 0, mol/L cell; 1, mol/L water; 2 mol/L rock
chem.SetUnitsPPassemblage(1);
// 0, mol/L cell; 1, mol/L water; 2 mol/L rock
chem.SetUnitsExchange(1);
// 0, mol/L cell; 1, mol/L water; 2 mol/L rock
chem.SetUnitsSurface(1);
// 0, mol/L cell; 1, mol/L water; 2 mol/L rock
chem.SetUnitsGasPhase(1);
// 0, mol/L cell; 1, mol/L water; 2 mol/L rock
chem.SetUnitsSSassemblage(1);
// 0, mol/L cell; 1, mol/L water; 2 mol/L rock
chem.SetUnitsKinetics(1);
// Set representative volume
std::vector<double> rv;
rv.resize(wp_size, 1.0);
chem.SetRepresentativeVolume(rv);
// Set initial porosity
std::vector<double> por;
por.resize(wp_size, 1);
chem.SetPorosity(por);
// Set initial saturation
std::vector<double> sat;
sat.resize(wp_size, 1.0);
chem.SetSaturation(sat);
// Load database
chem.LoadDatabase(database_path);
}
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 /* Iteration Count is dynamic, retrieving value from R (is only needed by
* master for the following loop) */ * master for the following loop) */
uint32_t maxiter = R.parseEval("mysetup$iterations"); uint32_t maxiter = params.timesteps.size();
double sim_time = .0; if (params.print_progressbar) {
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) {
chem.setProgressBarPrintout(true); chem.setProgressBarPrintout(true);
} }
@ -143,47 +259,33 @@ inline double RunMasterLoop(SimParams &params, RInside &R,
double dSimTime{0}; double dSimTime{0};
for (uint32_t iter = 1; iter < maxiter + 1; iter++) { for (uint32_t iter = 1; iter < maxiter + 1; iter++) {
double start_t = MPI_Wtime(); 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>( const double &dt = params.timesteps[iter - 1];
R.parseEval("mysetup$timesteps[" + std::to_string(iter) + "]"));
// cout << "CPP: Next time step is " << dt << "[s]" << endl; // cout << "CPP: Next time step is " << dt << "[s]" << endl;
MSG("Next time step is " + std::to_string(dt) + " [s]"); MSG("Next time step is " + std::to_string(dt) + " [s]");
/* displaying iteration number, with C++ and R iterator */ /* displaying iteration number, with C++ and R iterator */
MSG("Going through iteration " + std::to_string(iter)); MSG("Going through iteration " + std::to_string(iter));
MSG("R's $iter: " +
std::to_string((uint32_t)(R.parseEval("mysetup$iter"))) +
". Iteration");
/* run transport */ /* run transport */
// TODO: transport to diffusion
diffusion.simulate(dt); diffusion.simulate(dt);
chem.getField().update(diffusion.getField()); chem.getField().update(diffusion.getField());
// chem.getfield().update(diffusion.getfield());
MSG("Chemistry step"); MSG("Chemistry step");
chem.SetTimeStep(dt); chem.simulate(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");
// MDL master_iteration_end just writes on disk state_T and // MDL master_iteration_end just writes on disk state_T and
// state_C after every iteration if the cmdline option // state_C after every iteration if the cmdline option
// --ignore-results is not given (and thus the R variable // --ignore-results is not given (and thus the R variable
// store_result is TRUE) // 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) + "/" + MSG("End of *coupling* iteration " + std::to_string(iter) + "/" +
std::to_string(maxiter)); std::to_string(maxiter));
@ -194,175 +296,146 @@ inline double RunMasterLoop(SimParams &params, RInside &R,
dSimTime += end_t - start_t; dSimTime += end_t - start_t;
} // END SIMULATION LOOP } // 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(); Rcpp::List diffusion_profiling;
R.parseEvalQ("profiling$simtime_chemistry <- simtime_chemistry"); diffusion_profiling["simtime"] = diffusion.getTransportTime();
R["chemistry_loop"] = chem.GetMasterLoopTime(); if (params.use_dht) {
R.parseEvalQ("profiling$chemistry_loop <- chemistry_loop"); chem_profiling["dht_hits"] = Rcpp::wrap(chem.GetWorkerDHTHits());
chem_profiling["dht_evictions"] = Rcpp::wrap(chem.GetWorkerDHTEvictions());
R["chemistry_sequential"] = chem.GetMasterSequentialTime(); chem_profiling["dht_get_time"] = Rcpp::wrap(chem.GetWorkerDHTGetTimings());
R.parseEvalQ("profiling$simtime_sequential <- chemistry_sequential"); chem_profiling["dht_fill_time"] =
Rcpp::wrap(chem.GetWorkerDHTFillTimings());
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.getChemParams().use_interp) {
R["interp_w"] = Rcpp::wrap(chem.GetWorkerInterpolationWriteTimings()); if (params.use_interp) {
R.parseEvalQ("profiling$interp_write <- interp_w"); chem_profiling["interp_w"] =
R["interp_r"] = Rcpp::wrap(chem.GetWorkerInterpolationReadTimings()); Rcpp::wrap(chem.GetWorkerInterpolationWriteTimings());
R.parseEvalQ("profiling$interp_read <- interp_r"); chem_profiling["interp_r"] =
R["interp_g"] = Rcpp::wrap(chem.GetWorkerInterpolationGatherTimings()); Rcpp::wrap(chem.GetWorkerInterpolationReadTimings());
R.parseEvalQ("profiling$interp_gather <- interp_g"); chem_profiling["interp_g"] =
R["interp_fc"] = Rcpp::wrap(chem.GetWorkerInterpolationGatherTimings());
chem_profiling["interp_fc"] =
Rcpp::wrap(chem.GetWorkerInterpolationFunctionCallTimings()); Rcpp::wrap(chem.GetWorkerInterpolationFunctionCallTimings());
R.parseEvalQ("profiling$interp_function_calls <- interp_fc"); chem_profiling["interp_calls"] =
R["interp_calls"] = Rcpp::wrap(chem.GetWorkerInterpolationCalls()); Rcpp::wrap(chem.GetWorkerInterpolationCalls());
R.parseEvalQ("profiling$interp_calls <- interp_calls"); chem_profiling["interp_cached"] = Rcpp::wrap(chem.GetWorkerPHTCacheHits());
R["interp_cached"] = Rcpp::wrap(chem.GetWorkerPHTCacheHits());
R.parseEvalQ("profiling$interp_cached <- interp_cached");
} }
Rcpp::List profiling;
profiling["simtime"] = dSimTime;
profiling["chemistry"] = chem_profiling;
profiling["diffusion"] = diffusion_profiling;
chem.MasterLoopBreak(); chem.MasterLoopBreak();
diffusion.end();
return dSimTime; return profiling;
} }
int main(int argc, char *argv[]) { int main(int argc, char *argv[]) {
double dSimTime, sim_end; int world_size;
int world_size, world_rank;
MPI_Init(&argc, &argv); 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 (MY_RANK == 0) {
MSG("Running POET version " + std::string(poet_version));
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();
} }
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); 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*/ MSG("finished, cleanup of process " + std::to_string(MY_RANK));
// TODO: kann raus
R.parseEvalQ(kin_r_library);
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(); MPI_Finalize();
if (world_rank == 0) { if (MY_RANK == 0) {
MSG("done, bye!"); 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 <string>
#include <vector>
#include <Rcpp.h>
static const char *poet_version = "@POET_VERSION@"; static const char *poet_version = "@POET_VERSION@";
// using the Raw string literal to avoid escaping the quotes // 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 CONFIGURE_DEPENDS
"*.cpp" "*.c") "*.cpp" "*.c")
find_package(doctest REQUIRED)
add_executable(testPOET ${test_SRC}) 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") target_include_directories(testPOET PRIVATE "${PROJECT_SOURCE_DIR}/src")
get_filename_component(TEST_RInsideSourceFile "RInsidePOET_funcs.R" REALPATH) 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) { simple_field <- function(field) {
field <- as.data.frame(field, check.names = FALSE)
field$Na <- 0 field$Na <- 0
return(field) return(field)
} }
extended_field <- function(field, additional) { extended_field <- function(field, additional) {
field <- as.data.frame(field, check.names = FALSE)
additional <- as.data.frame(additional, check.names = FALSE)
return(field + additional) return(field + additional)
} }

5
test/testField.cpp
View File

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

3
test/testNamedVector.cpp
View File

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