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 'dht-rounding' into 'main'

Browse Source 
Merge latest changes including DHT rounding and code refactoring

See merge request sec34/port!14
This commit is contained in:
Max Lübke 2023-01-24 16:11:37 +01:00
commit 28bf539bec
24 changed files with 1318 additions and 265 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists.txt
View File

@ -22,6 +22,7 @@ add_subdirectory(src)
add_subdirectory(R_lib)
add_subdirectory(data)
add_subdirectory(app)
add_subdirectory(bench/dolo_diffu_inner)
add_subdirectory(ext/tug EXCLUDE_FROM_ALL)
add_subdirectory(ext/phreeqcrm EXCLUDE_FROM_ALL)

112
README.md
View File

@ -1,5 +1,5 @@
<!--
Time-stamp: "Last modified 2021-02-08 13:46:00 mluebke"
Time-stamp: "Last modified 2023-01-19 12:06:10 delucia"
-->
**Po**tsdamer **R**eactive **T**ransport
@ -7,7 +7,12 @@
# Forked Project
PORT is a fork of [POET](https://doi.org/10.5281/zenodo.4757913) integrating a standalone component for transport computations and leveraging PHREEQC_RM as geochemical solver. The following README is also applicable for this project.
*PORT* is a fork of [POET](https://doi.org/10.5281/zenodo.4757913)
integrating a standalone component for transport computations and
leveraging PHREEQC_RM as geochemical solver. The following README is
also applicable for this project.
![PORT's Coupling Scheme](./docs/20221216_Scheme_PORT_en.svg)
# POET
@ -31,9 +36,11 @@ To compile POET you need several software to be installed:
- R language and environment
- CMake 3.9+
If you want to build documentation during compilation, `doxygen`and `graphviz` must be provided too.
If you want to build documentation during compilation, `doxygen`and
`graphviz` must be provided too.
The following R libraries must then be installed, which will get the needed dependencies automatically:
The following R libraries must then be installed, which will get the
needed dependencies automatically:
- [devtools](https://www.r-project.org/nosvn/pandoc/devtools.html)
- [Rcpp](https://cran.r-project.org/web/packages/Rcpp/index.html)
@ -43,17 +50,21 @@ The following R libraries must then be installed, which will get the needed depe
### Compiling source code
The generation of makefiles is done with CMake. If you obtained POET from git, you should be able to generate Makefiles by running
The generation of makefiles is done with CMake. If you obtained POET
from git, you should be able to generate Makefiles by running
```sh
mkdir build && cd build
cmake ..
```
This will create the directory `build` and processes the CMake files and generate Makefiles from it. You're now able to run `make` to start build
process.
This will create the directory `build` and processes the CMake files
and generate Makefiles from it. You're now able to run `make` to start
build process.
If POET was obtained from the official SVN repository or the redmine at <https://redmine.cs.uni-potsdam.de/projects/poet> the branch or tag to be used have to be set via
If POET was obtained from the official SVN repository or the redmine
at <https://redmine.cs.uni-potsdam.de/projects/poet> the branch or tag
to be used have to be set via
```sh
mkdir build && cd build
@ -64,20 +75,26 @@ where currently available branches/tags are:
- dev
If everything went well you'll find the executable at `build/src/poet`, but it is recommended to install the POET project structure to a desired `CMAKE_INSTALL_PREFIX` with `make install`.
If everything went well you'll find the executable at
`build/src/poet`, but it is recommended to install the POET project
structure to a desired `CMAKE_INSTALL_PREFIX` with `make install`.
During the generation of Makefiles, various options can be specified via `cmake -D <option>=<value> [...]`. Currently there are the following available options:
During the generation of Makefiles, various options can be specified
via `cmake -D <option>=<value> [...]`. Currently there are the
following available options:
- **DHT_Debug**=_boolean_ - toggles the output of detailed statistics about DHT
usage (`cmake -D DHT_Debug=ON`). Defaults to _OFF_.
- **BUILD_DOC**=_boolean_ - toggles the generation of documantiation during
compilation process. Defaults to _ON_.
- _only from svn version:_ **POET_SET_BRANCH**=_string_ - set branch or tag whose code is used
- **DHT_Debug**=_boolean_ - toggles the output of detailed statistics
about DHT usage (`cmake -D DHT_Debug=ON`). Defaults to _OFF_.
- **BUILD_DOC**=_boolean_ - toggles the generation of documantiation
during compilation process. Defaults to _ON_.
- _only from svn version:_ **POET_SET_BRANCH**=_string_ - set branch
or tag whose code is used
### Example: Build from scratch
Assuming that only the C/C++ compiler, MPI libraries, R runtime environment and
CMake have been installed, POET can be installed as follows:
Assuming that only the C/C++ compiler, MPI libraries, R runtime
environment and CMake have been installed, POET can be installed as
follows:
```sh
# start R environment
@ -99,9 +116,10 @@ $ make -j<max_numprocs>
$ make install
```
This will install a POET project structure into `/home/<user>/poet` which is
called hereinafter `<POET_INSTALL_DIR>`. With this version of POET we **do not
recommend** to install to hierarchies like `/usr/local/` etc.
This will install a POET project structure into `/home/<user>/poet`
which is called hereinafter `<POET_INSTALL_DIR>`. With this version of
POET we **do not recommend** to install to hierarchies like
`/usr/local/` etc.
The correspondending directory tree would look like this:
@ -122,21 +140,23 @@ The correspondending directory tree would look like this:
```
The R libraries will be loaded at runtime and the paths are hardcoded
absolute paths inside `poet.cpp`. So, if you consider to move `bin/poet` either
change paths of the R source files and recompile POET or also move `R_lib/*`
according to the binary.
absolute paths inside `poet.cpp`. So, if you consider to move
`bin/poet` either change paths of the R source files and recompile
POET or also move `R_lib/*` according to the binary.
To display the generated html documentation just open `docs/html/index.html`
with the browser of your choice.
To display the generated html documentation just open
`docs/html/index.html` with the browser of your choice.
## Running
Before POET is ready to run, a working directory must be created. In this
directory you should find the executable file, the R scripts
`<POET_ROOT>/R_lib/kin_r_library.R` and `<POET_ROOT>/R_lib/parallel_r_library.R`
and the simulation description e.g. `<POET_ROOT>/data/chem_problems/SimDol2D.R`.
Before POET is ready to run, a working directory must be created. In
this directory you should find the executable file, the R scripts
`<POET_ROOT>/R_lib/kin_r_library.R` and
`<POET_ROOT>/R_lib/parallel_r_library.R` and the simulation
description e.g. `<POET_ROOT>/data/chem_problems/SimDol2D.R`.
Run POET by `mpirun ./poet <OPTIONS> <SIMFILE> <OUTPUT_DIRECTORY>` where:
Run POET by `mpirun ./poet <OPTIONS> <SIMFILE> <OUTPUT_DIRECTORY>`
where:
- **OPTIONS** - runtime parameters (explained below)
- **SIMFILE** - simulation described as R script (currently supported:
@ -161,8 +181,8 @@ The following parameters can be set:
#### Additions to `dht-signif`
Only used if no vector is given in setup file. For individual values per column
use R vector `signif_vector` in `SIMFILE`.
Only used if no vector is given in setup file. For individual values
per column use R vector `signif_vector` in `SIMFILE`.
#### Additions to `dht-snaps`
@ -176,22 +196,24 @@ Following values can be set:
### Example: Running from scratch
We will continue the above example and start a simulation with `SimDol2D.R`,
which is the only simulation supported at this moment. The required flow velocities
snapshots are included in the R package Rmufits. It's a 2D, 50x50 grid, with 20 time
steps. To start the simulation with 4 processes `cd` into your previously installed
We will continue the above example and start a simulation with
`SimDol2D.R`, which is the only simulation supported at this moment.
The required flow velocities snapshots are included in the R package
Rmufits. It's a 2D, 50x50 grid, with 20 time steps. To start the
simulation with 4 processes `cd` into your previously installed
POET-dir `<POET_INSTALL_DIR>/bin` and run:
```sh
mpirun -n 4 ./poet ../data/SimDol2D.R output
```
After a finished simulation all data generated by POET will be found in the
directory `output`.
After a finished simulation all data generated by POET will be found
in the directory `output`.
You might want to use the DHT to cache previously simulated data and
reuse them in further time-steps. Just append `--dht` to the options of POET to
activate the usage of the DHT. The resulting call would look like this:
reuse them in further time-steps. Just append `--dht` to the options
of POET to activate the usage of the DHT. The resulting call would
look like this:
```sh
mpirun -n 4 ./poet --dht SimDol2D.R output
@ -199,9 +221,9 @@ mpirun -n 4 ./poet --dht SimDol2D.R output
## About the usage of MPI_Wtime()
Implemented time measurement functions uses `MPI_Wtime()`. Some important
information from the OpenMPI Man Page:
Implemented time measurement functions uses `MPI_Wtime()`. Some
important information from the OpenMPI Man Page:
For example, on platforms that support it, the clock_gettime() function will be
used to obtain a monotonic clock value with whatever precision is supported on
that platform (e.g., nanoseconds).
For example, on platforms that support it, the clock_gettime()
function will be used to obtain a monotonic clock value with whatever
precision is supported on that platform (e.g., nanoseconds).

48
R_lib/kin_r_library.R
View File

@ -46,6 +46,17 @@ master_init <- function(setup) {
setup$maxiter <- setup$iterations
setup$timesteps <- setup$timesteps
setup$simulation_time <- 0
if (!(exists("setup$store_result"))) {
setup$store_result <- TRUE
}
if (setup$store_result) {
if (!(exists("setup$out_save"))) {
setup$out_save <- seq(1, setup$iterations)
}
}
return(setup)
}
@ -56,22 +67,22 @@ master_iteration_end <- function(setup) {
iter <- setup$iter
## MDL Write on disk state_T and state_C after every iteration
## comprised in setup$out_save
# if (store_result) {
# if (iter %in% setup$out_save) {
nameout <- paste0(fileout, "/iter_", sprintf("%03d", iter), ".rds")
info <- list(
tr_req_dt = as.integer(setup$req_dt)
# tr_allow_dt = setup$allowed_dt,
# tr_inniter = as.integer(setup$inniter)
)
saveRDS(list(
T = setup$state_T, C = setup$state_C,
simtime = as.integer(setup$simtime),
tr_info = info
), file = nameout)
msgm("results stored in <", nameout, ">")
# }
# }
if (setup$store_result) {
if (iter %in% setup$out_save) {
nameout <- paste0(fileout, "/iter_", sprintf("%03d", iter), ".rds")
info <- list(
tr_req_dt = as.integer(setup$req_dt)
# tr_allow_dt = setup$allowed_dt,
# tr_inniter = as.integer(setup$inniter)
)
saveRDS(list(
T = setup$state_T, C = setup$state_C,
simtime = as.integer(setup$simtime),
tr_info = info
), file = nameout)
msgm("results stored in <", nameout, ">")
}
}
msgm("done iteration", iter, "/", setup$maxiter)
setup$iter <- setup$iter + 1
return(setup)
@ -264,3 +275,8 @@ StoreSetup <- function(setup) {
saveRDS(to_store, file = paste0(fileout, "/setup.rds"))
msgm("initialization stored in ", paste0(fileout, "/setup.rds"))
}
GetWorkPackageSizesVector <- function(n_packages, package_size, len) {
ids <- rep(1:n_packages, times=package_size, each = 1)[1:len]
return(as.integer(table(ids)))
}

6
bench/dolo_diffu_inner/CMakeLists.txt Normal file
View File

@ -0,0 +1,6 @@
install(FILES
dolo_diffu_inner.R
dolo_inner.pqi
DESTINATION
share/poet/bench
)

51
bench/dolo_diffu_inner/Eval.R Normal file
View File

@ -0,0 +1,51 @@
## 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()

146
bench/dolo_diffu_inner/dolo_diffu_inner.R Normal file
View File

@ -0,0 +1,146 @@
## Time-stamp: "Last modified 2023-01-10 13:51:40 delucia"
database <- normalizePath("../share/poet/examples/phreeqc_kin.dat")
input_script <- normalizePath("../share/poet/bench/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" = 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],
init_cell = as.data.frame(init_cell),
props = names(init_cell),
database = database,
input_script = input_script
)
##################################################################
## Section 2 ##
## Diffusion parameters and boundary conditions ##
##################################################################
## initial conditions
init_diffu <- c(
"H" = 110.683,
"O" = 55.3413,
"Charge" = -5.0822e-19,
"C" = 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" = 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" = 0,
"Ca" = 0,
"Cl" = 0.002,
"Mg" = 0.001
),
list(
"H" = 110.683,
"O" = 55.3413,
"Charge" = 1.90431e-16,
"C" = 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)
diffusion <- list(
init = init_diffu,
vecinj = do.call(rbind.data.frame, vecinj_diffu),
vecinj_inner = vecinj_inner,
vecinj_index = boundary,
alpha = alpha_diffu
)
#################################################################
## Section 3 ##
## Chemistry module (Phreeqc) ##
#################################################################
## # Needed when using DHT
signif_vector <- c(10, 10, 2, 5, 5, 5, 5, 0, 5, 5)
prop_type <- c("", "", "", "act", "act", "act", "act", "ignore", "", "")
prop <- names(init_cell)
chemistry <- list(
database = database,
input_script = input_script
)
#################################################################
## Section 4 ##
## Putting all those things together ##
#################################################################
iterations <- 1000
dt <- 200
setup <- list(
grid = grid,
diffusion = diffusion,
chemistry = chemistry,
iterations = iterations,
timesteps = rep(dt, iterations),
store_result = TRUE,
out_save = c(5, iterations)
)

146
bench/dolo_diffu_inner/dolo_diffu_inner_large.R Normal file
View File

@ -0,0 +1,146 @@
## Time-stamp: "Last modified 2023-01-10 13:51:40 delucia"
database <- normalizePath("../share/poet/examples/phreeqc_kin.dat")
input_script <- normalizePath("../share/poet/bench/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(2, 1),
type = types[1],
init_cell = as.data.frame(init_cell),
props = names(init_cell),
database = database,
input_script = input_script
)
##################################################################
## Section 2 ##
## Diffusion parameters and boundary conditions ##
##################################################################
## initial conditions
init_diffu <- c(
"H" = 110.683,
"O" = 55.3413,
"Charge" = -5.0822e-19,
"C" = 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" = 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" = 0,
"Ca" = 0,
"Cl" = 0.002,
"Mg" = 0.001
),
list(
"H" = 110.683,
"O" = 55.3413,
"Charge" = 1.90431e-16,
"C" = 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)
diffusion <- list(
init = init_diffu,
vecinj = do.call(rbind.data.frame, vecinj_diffu),
vecinj_inner = vecinj_inner,
vecinj_index = boundary,
alpha = alpha_diffu
)
#################################################################
## Section 3 ##
## Chemistry module (Phreeqc) ##
#################################################################
## # Needed when using DHT
signif_vector <- c(10, 10, 2, 5, 5, 5, 5, 0, 5, 5)
prop_type <- c("", "", "", "act", "act", "act", "act", "ignore", "", "")
prop <- names(init_cell)
chemistry <- list(
database = database,
input_script = input_script
)
#################################################################
## 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 = c(5, iterations)
)

35
bench/dolo_diffu_inner/dolo_inner.pqi Normal file
View File

@ -0,0 +1,35 @@
SELECTED_OUTPUT
-high_precision true
-reset false
-time
-soln
-temperature true
-water true
-pH
-pe
-totals C Ca Cl Mg
-kinetic_reactants Calcite Dolomite
-equilibrium O2g
SOLUTION 1
units mol/kgw
temp 25.0
water 1
pH 9.91 charge
pe 4.0
C 1.2279E-04
Ca 1.2279E-04
Cl 1E-12
Mg 1E-12
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

2
data/CMakeLists.txt
View File

@ -5,4 +5,4 @@ FILES
phreeqc_kin.dat
dol.pqi
DESTINATION
data)
share/poet/examples)

12
data/SimDol2D_diffu.R
View File

@ -1,13 +1,13 @@
database <- normalizePath("../data/phreeqc_kin.dat")
input_script <- normalizePath("../data/dol.pqi")
database <- normalizePath("../share/poet/examples/phreeqc_kin.dat")
input_script <- normalizePath("../share/poet/examples/dol.pqi")
#################################################################
## Section 1 ##
## Grid initialization ##
#################################################################
n <- 5
m <- 5
n <- 100
m <- 100
types <- c("scratch", "phreeqc", "rds")
@ -168,8 +168,8 @@ selout <- c(
# Needed when using DHT
signif_vector <- c(7, 7, 7, 7, 7, 7, 7, 5, 5, 7)
prop_type <- c("act", "act", "act", "act", "logact", "logact", "ignore", "act", "act", "act")
signif_vector <- c(10, 10, 10, 7, 7, 7, 7, 0, 5, 5)
prop_type <- c("", "", "", "act", "act", "act", "act", "ignore", "", "")
prop <- names(init_cell)
chemistry <- list(

572
docs/20221216_Scheme_PORT_en.svg Normal file
View File

File diff suppressed because one or more lines are too long
Side by Side

After

Width:  |  Height:  |  Size: 73 KiB

2
ext/tug

@ -1 +1 @@
Subproject commit d4e3ab8544628c72921bf60d4462d6d131b0ffb5
Subproject commit 79d7a32fc25fd97b78f320bba854b1db64e059f7

47
include/poet/DHT_Wrapper.hpp
View File

@ -33,16 +33,23 @@ extern "C" {
#include <mpi.h>
/**
* @brief Cut double value after signif digit
*
* Macro to round a double value by cutting every digit after significant digit
*
*/
// #define ROUND(value, signif) \
// (((int)(pow(10.0, (double)signif) * value)) * pow(10.0, (double)-signif))
namespace poet {
using DHT_Keyelement = struct keyelem {
std::int8_t exp : 8;
std::int64_t significant : 56;
};
using DHT_ResultObject = struct DHTResobj {
uint32_t length;
std::vector<std::vector<DHT_Keyelement>> keys;
std::vector<std::vector<double>> results;
std::vector<bool> needPhreeqc;
void ResultsToMapping(std::vector<int32_t> &curr_mapping);
void ResultsToWP(std::vector<double> &curr_wp);
};
/**
* @brief Return user-defined md5sum
*
@ -76,11 +83,11 @@ public:
* @param dht_comm Communicator which addresses all participating DHT
* processes
* @param buckets_per_process Count of buckets to allocate for each process
* @param data_size Size of data in bytes
* @param key_size Size of key in bytes
* @param key_count Count of key entries
* @param data_count Count of data entries
*/
DHT_Wrapper(SimParams &params, MPI_Comm dht_comm, int buckets_per_process,
int data_size, int key_size);
DHT_Wrapper(const poet::SimParams &params, MPI_Comm dht_comm,
uint32_t dht_size, uint32_t key_count, uint32_t data_count);
/**
* @brief Destroy the dht wrapper object
*
@ -108,9 +115,8 @@ public:
* @param[in,out] work_package Pointer to current work package
* @param dt Current timestep of simulation
*/
auto checkDHT(int length, std::vector<bool> &out_result_index,
double *work_package, double dt)
-> std::vector<std::vector<double>>;
auto checkDHT(int length, double dt, const std::vector<double> &work_package)
-> poet::DHT_ResultObject;
/**
* @brief Write simulated values into DHT
@ -128,8 +134,8 @@ public:
* outputs of the PHREEQC simulation
* @param dt Current timestep of simulation
*/
void fillDHT(int length, std::vector<bool> &result_index,
double *work_package, double *results, double dt);
void fillDHT(int length, const DHT_ResultObject &dht_check_data,
const std::vector<double> &results);
/**
* @brief Dump current DHT state into file.
@ -184,6 +190,9 @@ public:
uint64_t getEvictions();
private:
uint32_t key_count;
uint32_t data_count;
/**
* @brief Transform given workpackage into DHT key
*
@ -211,7 +220,7 @@ private:
* @param key Pointer to work package handled as the key
* @param dt Current time step of the simulation
*/
std::vector<double> fuzzForDHT(int var_count, void *key, double dt);
std::vector<DHT_Keyelement> fuzzForDHT(int var_count, void *key, double dt);
/**
* @brief DHT handle

3
include/poet/DiffusionModule.hpp
View File

@ -24,6 +24,7 @@
#include "poet/SimParams.hpp"
#include <array>
#include <bits/stdint-uintn.h>
#include <cmath>
#include <poet/Grid.hpp>
#include <string>
#include <tug/BoundaryCondition.hpp>
@ -85,6 +86,8 @@ private:
void initialize(poet::DiffusionParams args);
void RoundToZero(double *field, uint32_t cell_count) const;
Grid &grid;
uint8_t dim;

13
include/poet/SimParams.hpp
View File

@ -23,6 +23,7 @@
#include <bits/stdint-uintn.h>
#include <string>
#include <string_view>
#include <vector>
#include "argh.hpp" // Argument handler https://github.com/adishavit/argh
@ -184,7 +185,7 @@ public:
*
* @return t_simparams Parameter struct
*/
t_simparams getNumParams();
t_simparams getNumParams() const;
/**
* @brief Get the DHT_Signif_Vector
@ -195,7 +196,7 @@ public:
* @return std::vector<int> Vector of integers containing information about
* significant digits
*/
std::vector<int> getDHTSignifVector();
std::vector<int> getDHTSignifVector() const;
/**
* @brief Get the DHT_Prop_Type_Vector
@ -205,7 +206,7 @@ public:
* @return std::vector<std::string> Vector if strings defining a type of a
* variable
*/
std::vector<std::string> getDHTPropTypeVector();
std::vector<std::string> getDHTPropTypeVector() const;
/**
* @brief Return name of DHT snapshot.
@ -215,7 +216,7 @@ public:
*
* @return std::string Absolute paht to the DHT snapshot
*/
std::string getDHTFile();
std::string_view getDHTFile() const;
/**
* @brief Get the filesim name
@ -225,7 +226,7 @@ public:
*
* @return std::string Absolute path to R file
*/
std::string getFilesim();
std::string_view getFilesim() const;
/**
* @brief Get the output directory
@ -235,7 +236,7 @@ public:
*
* @return std::string Absolute path to output path
*/
std::string getOutDir();
std::string_view getOutDir() const;
private:
/**

6
src/CMakeLists.txt
View File

@ -5,8 +5,14 @@ file(GLOB poet_lib_SRC
find_library(MATH_LIBRARY m)
find_library(CRYPTO_LIBRARY crypto)
option(POET_DHT_DEBUG "Build with DHT debug info" OFF)
add_library(poet_lib ${poet_lib_SRC})
target_include_directories(poet_lib PUBLIC ${PROJECT_SOURCE_DIR}/include)
target_link_libraries(poet_lib PUBLIC
MPI::MPI_C ${MATH_LIBRARY} ${CRYPTO_LIBRARY} RRuntime tug PhreeqcRM)
target_compile_definitions(poet_lib PUBLIC STRICT_R_HEADERS OMPI_SKIP_MPICXX)
if(POET_DHT_DEBUG)
target_compile_definitions(poet_lib PRIVATE DHT_STATISTICS)
endif()

11
src/ChemMaster.cpp
View File

@ -26,6 +26,7 @@
#include <Rcpp.h>
#include <array>
#include <cstdint>
#include <cstring>
#include <iostream>
#include <string>
@ -56,13 +57,9 @@ ChemMaster::ChemMaster(SimParams &params, RInside &R_, Grid &grid_)
uint32_t n_packages = (uint32_t)(grid.GetTotalCellCount() / this->wp_size) +
(mod_pkgs != 0 ? 1 : 0);
this->wp_sizes_vector = std::vector<uint32_t>(n_packages, this->wp_size);
if (mod_pkgs) {
auto itEndVector = this->wp_sizes_vector.end() - 1;
for (uint32_t i = 0; i < this->wp_size - mod_pkgs; i++) {
*(itEndVector - i) -= 1;
}
}
Rcpp::Function wp_f("GetWorkPackageSizesVector");
this->wp_sizes_vector = Rcpp::as<std::vector<uint32_t>>(
wp_f(n_packages, this->wp_size, grid.GetTotalCellCount()));
this->state = this->grid.RegisterState(
poet::BaseChemModule::CHEMISTRY_MODULE_NAME, this->prop_names);

53
src/ChemWorker.cpp
View File

@ -19,10 +19,12 @@
*/
#include "poet/ChemSimPar.hpp"
#include "poet/DHT_Wrapper.hpp"
#include "poet/SimParams.hpp"
#include <Rcpp.h>
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <mpi.h>
@ -57,20 +59,20 @@ ChemWorker::ChemWorker(SimParams &params, RInside &R_, Grid &grid_,
<< endl;
if (this->dht_enabled) {
int data_size = this->prop_names.size() * sizeof(double);
int key_size = (this->prop_names.size() - 1) * sizeof(double) +
(dt_differ * sizeof(double));
int dht_buckets_per_process =
dht_size_per_process / (1 + data_size + key_size);
uint32_t iKeyCount = this->prop_names.size() + (dt_differ);
uint32_t iDataCount = this->prop_names.size();
if (world_rank == 1)
cout << "CPP: Worker: data size: " << data_size << " bytes" << endl
<< "CPP: Worker: key size: " << key_size << " bytes" << endl
<< "CPP: Worker: buckets per process " << dht_buckets_per_process
<< endl;
cout << "CPP: Worker: data count: " << iDataCount << " entries" << endl
<< "CPP: Worker: key count: " << iKeyCount << " entries" << endl
<< "CPP: Worker: memory per process "
<< params.getNumParams().dht_size_per_process / std::pow(10, 6)
<< " MByte" << endl;
dht = new DHT_Wrapper(params, dht_comm, dht_buckets_per_process, data_size,
key_size);
dht = new DHT_Wrapper(params, dht_comm,
params.getNumParams().dht_size_per_process, iKeyCount,
iDataCount);
if (world_rank == 1)
cout << "CPP: Worker: DHT created!" << endl;
@ -171,11 +173,9 @@ void ChemWorker::doWork(MPI_Status &probe_status) {
std::vector<double> vecCurrWP(
mpi_buffer,
mpi_buffer + (local_work_package_size * this->prop_names.size()));
std::vector<int32_t> vecMappingWP(this->wp_size);
std::vector<std::vector<double>> vecDHTResults;
std::vector<double> vecDHTKeys;
std::vector<bool> vecNeedPhreeqc(this->wp_size, true);
DHT_ResultObject DHT_Results;
for (uint32_t i = 0; i < local_work_package_size; i++) {
vecMappingWP[i] = i;
@ -195,19 +195,10 @@ void ChemWorker::doWork(MPI_Status &probe_status) {
if (dht_enabled) {
/* check for values in DHT */
dht_get_start = MPI_Wtime();
vecDHTKeys = this->phreeqc_rm->ReplaceTotalsByPotentials(
vecCurrWP, local_work_package_size);
vecDHTResults = dht->checkDHT(local_work_package_size, vecNeedPhreeqc,
vecDHTKeys.data(), dt);
DHT_Results = dht->checkDHT(local_work_package_size, dt, vecCurrWP);
dht_get_end = MPI_Wtime();
uint32_t iMappingIndex = 0;
for (uint32_t i = 0; i < local_work_package_size; i++) {
if (vecMappingWP[i] == -1) {
continue;
}
vecMappingWP[i] = (vecNeedPhreeqc[i] ? iMappingIndex++ : -1);
}
DHT_Results.ResultsToMapping(vecMappingWP);
}
phreeqc_time_start = MPI_Wtime();
@ -216,12 +207,7 @@ void ChemWorker::doWork(MPI_Status &probe_status) {
phreeqc_time_end = MPI_Wtime();
if (dht_enabled) {
for (uint32_t i = 0; i < local_work_package_size; i++) {
if (!vecNeedPhreeqc[i]) {
std::copy(vecDHTResults[i].begin(), vecDHTResults[i].end(),
vecCurrWP.begin() + (this->prop_names.size() * i));
}
}
DHT_Results.ResultsToWP(vecCurrWP);
}
/* send results to master */
@ -232,8 +218,7 @@ void ChemWorker::doWork(MPI_Status &probe_status) {
if (dht_enabled) {
/* write results to DHT */
dht_fill_start = MPI_Wtime();
dht->fillDHT(local_work_package_size, vecNeedPhreeqc, vecDHTKeys.data(),
vecCurrWP.data(), dt);
dht->fillDHT(local_work_package_size, DHT_Results, vecCurrWP);
dht_fill_end = MPI_Wtime();
timing[1] += dht_get_end - dht_get_start;

194
src/DHT_Wrapper.cpp
View File

@ -19,9 +19,11 @@
*/
#include "poet/HashFunctions.hpp"
#include <cstddef>
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <cstring>
#include <openssl/evp.h>
#include <poet/DHT_Wrapper.hpp>
@ -33,19 +35,49 @@
using namespace poet;
using namespace std;
inline double round_signif(double value, int32_t signif) {
const double multiplier = std::pow(10.0, signif);
return .0 + std::trunc(value * multiplier) / multiplier;
inline DHT_Keyelement round_key_element(double value, std::uint32_t signif) {
std::int8_t exp = (std::int8_t)std::floor(std::log10(std::fabs(value)));
std::int64_t significant = value * std::pow(10, signif - exp);
DHT_Keyelement elem;
elem.exp = exp;
elem.significant = significant;
return elem;
}
DHT_Wrapper::DHT_Wrapper(SimParams &params, MPI_Comm dht_comm,
int buckets_per_process, int data_size, int key_size) {
void poet::DHT_ResultObject::ResultsToMapping(
std::vector<int32_t> &curr_mapping) {
uint32_t iMappingIndex = 0;
for (uint32_t i = 0; i < this->length; i++) {
if (curr_mapping[i] == -1) {
continue;
}
curr_mapping[i] = (this->needPhreeqc[i] ? iMappingIndex++ : -1);
}
}
void poet::DHT_ResultObject::ResultsToWP(std::vector<double> &curr_wp) {
for (uint32_t i = 0; i < this->length; i++) {
if (!this->needPhreeqc[i]) {
const uint32_t length = this->results[i].end() - this->results[i].begin();
std::copy(this->results[i].begin(), this->results[i].end(),
curr_wp.begin() + (length * i));
}
}
}
DHT_Wrapper::DHT_Wrapper(const poet::SimParams &params, MPI_Comm dht_comm,
uint32_t dht_size, uint32_t key_count,
uint32_t data_count)
: key_count(key_count), data_count(data_count) {
poet::initHashCtx(EVP_md5());
// initialize DHT object
uint32_t key_size = key_count * sizeof(DHT_Keyelement);
uint32_t data_size = data_count * sizeof(double);
uint32_t buckets_per_process = dht_size / (1 + data_size + key_size);
dht_object = DHT_create(dht_comm, buckets_per_process, data_size, key_size,
&poet::hashDHT);
// allocate memory for fuzzing buffer
fuzzing_buffer = (double *)malloc(key_size);
// extract needed values from sim_param struct
t_simparams tmp = params.getNumParams();
@ -64,76 +96,61 @@ DHT_Wrapper::~DHT_Wrapper() {
free(fuzzing_buffer);
poet::freeHashCtx();
}
auto DHT_Wrapper::checkDHT(int length, double dt,
const std::vector<double> &work_package)
-> poet::DHT_ResultObject {
DHT_ResultObject check_data;
check_data.length = length;
check_data.keys.resize(length);
check_data.results.resize(length);
check_data.needPhreeqc.resize(length);
auto DHT_Wrapper::checkDHT(int length, std::vector<bool> &out_result_index,
double *work_package, double dt)
-> std::vector<std::vector<double>> {
void *key;
int res;
// var count -> count of variables per grid cell
int var_count = dht_prop_type_vector.size() - 1;
std::vector<std::vector<double>> data(length);
// loop over every grid cell contained in work package
for (int i = 0; i < length; i++) {
// point to current grid cell
key = (void *)&(work_package[i * var_count]);
data[i].resize(dht_prop_type_vector.size());
void *key = (void *)&(work_package[i * this->key_count]);
auto &data = check_data.results[i];
auto &key_vector = check_data.keys[i];
// fuzz data (round, logarithm etc.)
std::vector<double> vecFuzz = fuzzForDHT(var_count, key, dt);
data.resize(this->data_count);
key_vector = fuzzForDHT(this->key_count, key, dt);
// overwrite input with data from DHT, IF value is found in DHT
res = DHT_read(dht_object, vecFuzz.data(), data[i].data());
int res = DHT_read(this->dht_object, key_vector.data(), data.data());
// if DHT_SUCCESS value was found ...
if (res == DHT_SUCCESS) {
// ... and grid cell will be marked as 'not to be simulating'
out_result_index[i] = false;
dht_hits++;
}
// ... otherwise ...
else if (res == DHT_READ_MISS) {
// grid cell needs to be simulated by PHREEQC
dht_miss++;
} else {
// MPI ERROR ... WHAT TO DO NOW?
// RUNNING CIRCLES WHILE SCREAMING
switch (res) {
case DHT_SUCCESS:
check_data.needPhreeqc[i] = false;
this->dht_hits++;
break;
case DHT_READ_MISS:
check_data.needPhreeqc[i] = true;
this->dht_miss++;
break;
}
}
return data;
return check_data;
}
void DHT_Wrapper::fillDHT(int length, std::vector<bool> &result_index,
double *work_package, double *results, double dt) {
void *key;
void *data;
int res;
// var count -> count of variables per grid cell
int var_count = dht_prop_type_vector.size();
void DHT_Wrapper::fillDHT(int length, const DHT_ResultObject &dht_check_data,
const std::vector<double> &results) {
// loop over every grid cell contained in work package
for (int i = 0; i < length; i++) {
key = (void *)&(work_package[i * (var_count - 1)]);
data = (void *)&(results[i * var_count]);
// If true grid cell was simulated, needs to be inserted into dht
if (result_index[i]) {
if (dht_check_data.needPhreeqc[i]) {
const auto &key = dht_check_data.keys[i];
void *data = (void *)&(results[i * this->data_count]);
// fuzz data (round, logarithm etc.)
std::vector<double> vecFuzz = fuzzForDHT(var_count - 1, key, dt);
// insert simulated data with fuzzed key into DHT
res = DHT_write(dht_object, vecFuzz.data(), data);
int res = DHT_write(this->dht_object, (void *)key.data(), data);
// if data was successfully written ...
if (res != DHT_SUCCESS) {
// ... also check if a previously written value was evicted
if (res == DHT_WRITE_SUCCESS_WITH_EVICTION) {
// and increment internal eviciton counter
dht_evictions++;
} else {
// MPI ERROR ... WHAT TO DO NOW?
// RUNNING CIRCLES WHILE SCREAMING
}
if ((res != DHT_SUCCESS) && (res == DHT_WRITE_SUCCESS_WITH_EVICTION)) {
dht_evictions++;
}
}
}
@ -173,67 +190,32 @@ uint64_t DHT_Wrapper::getMisses() { return this->dht_miss; }
uint64_t DHT_Wrapper::getEvictions() { return this->dht_evictions; }
std::vector<double> DHT_Wrapper::fuzzForDHT(int var_count, void *key,
double dt) {
std::vector<DHT_Keyelement> DHT_Wrapper::fuzzForDHT(int var_count, void *key,
double dt) {
constexpr double zero_val = 10E-14;
std::vector<DHT_Keyelement> vecFuzz(var_count);
std::memset(&vecFuzz[0], 0, sizeof(DHT_Keyelement) * var_count);
std::vector<double> vecFuzz(var_count, .0);
unsigned int i = 0;
// introduce fuzzing to allow more hits in DHT
// loop over every variable of grid cell
for (i = 0; i < (unsigned int)var_count; i++) {
// check if variable is defined as 'act'
if (dht_prop_type_vector[i] == "act") {
// if log is enabled (default)
if (dht_log) {
// if variable is smaller than 0, which would be a strange result,
// warn the user and set fuzzing_buffer to 0 at this index
if (((double *)key)[i] < 0) {
cerr << "dht_wrapper.cpp::fuzz_for_dht(): Warning! Negative value in "
"key!"
<< endl;
vecFuzz[i] = 0;
}
// if variable is 0 set fuzzing buffer to 0
else if (((double *)key)[i] == 0)
vecFuzz[i] = 0;
// otherwise ...
else {
// round current variable value by applying log with base 10, negate
// (since the actual values will be between 0 and 1) and cut result
// after significant digit
vecFuzz[i] = round_signif(-(std::log10(((double *)key)[i])),
dht_signif_vector[i]);
}
double &curr_key = ((double *)key)[i];
if (curr_key != 0) {
if (curr_key < zero_val && this->dht_prop_type_vector[i] == "act") {
continue;
}
// if log is disabled
else {
// just round by cutting after signifanct digit
vecFuzz[i] = round_signif((((double *)key)[i]), dht_signif_vector[i]);
if (this->dht_prop_type_vector[i] == "ignore") {
continue;
}
} else if (dht_prop_type_vector[i] == "charge") {
vecFuzz[i] = round_signif(((double *)key)[i], dht_signif_vector[i]);
}
// if variable is defined as 'logact' (log was already applied e.g. pH)
else if (dht_prop_type_vector[i] == "logact") {
// just round by cutting after signifanct digit
vecFuzz[i] = round_signif((((double *)key)[i]), dht_signif_vector[i]);
}
// if defined ass 'ignore' ...
else if (dht_prop_type_vector[i] == "ignore") {
// ... just set fuzzing buffer to 0
vecFuzz[i] = 0;
}
// and finally, if type is not defined, print error message
else {
cerr << "dht_wrapper.cpp::fuzz_for_dht(): Warning! Probably wrong "
"prop_type!"
<< endl;
vecFuzz[i] = round_key_element(curr_key, dht_signif_vector[i]);
}
}
// if timestep differs over iterations set current current time step at the
// end of fuzzing buffer
if (dt_differ) {
vecFuzz[var_count] = dt;
vecFuzz[var_count] = round_key_element(dt, 55);
}
return vecFuzz;

14
src/DiffusionModule.cpp
View File

@ -23,6 +23,7 @@
#include "tug/Diffusion.hpp"
#include <Rcpp.h>
#include <algorithm>
#include <bits/stdint-intn.h>
#include <cstdint>
#include <poet/ChemSimSeq.hpp>
#include <poet/DiffusionModule.hpp>
@ -36,6 +37,8 @@
using namespace poet;
static constexpr double ZERO_MULTIPLICATOR = 10E-14;
constexpr std::array<uint8_t, 4> borders = {
tug::bc::BC_SIDE_LEFT, tug::bc::BC_SIDE_RIGHT, tug::bc::BC_SIDE_TOP,
tug::bc::BC_SIDE_BOTTOM};
@ -180,6 +183,11 @@ void DiffusionModule::simulate(double dt) {
} else {
tug::diffusion::ADI_2D(this->diff_input, in_field, in_alpha.data());
}
// TODO: do not use hardcoded index for O, H and charge
if (i > 2) {
this->RoundToZero(in_field, this->n_cells_per_prop);
}
}
std::cout << " done!\n";
@ -188,6 +196,12 @@ void DiffusionModule::simulate(double dt) {
transport_t += sim_a_transport - sim_b_transport;
}
inline void DiffusionModule::RoundToZero(double *field,
uint32_t cell_count) const {
for (uint32_t i = 0; i < cell_count; i++) {
field[i] = ((int32_t)(field[i] / ZERO_MULTIPLICATOR)) * ZERO_MULTIPLICATOR;
}
}
void DiffusionModule::end() {
// R["simtime_transport"] = transport_t;

5
src/PhreeqcWrapper.cpp
View File

@ -52,7 +52,7 @@ void poet::PhreeqcWrapper::SetupAndLoadDB(
// Set initial porosity
std::vector<double> por;
por.resize(this->iWPSize, 0.05);
por.resize(this->iWPSize, 1);
this->SetPorosity(por);
// Set initial saturation
@ -66,7 +66,8 @@ void poet::PhreeqcWrapper::SetupAndLoadDB(
void poet::PhreeqcWrapper::InitFromFile(const std::string &strInputFile) {
this->RunFile(true, true, false, strInputFile);
this->RunString(true, false, true, "DELETE; -all");
// MDL: this is run only by the workers
this->RunString(true, false, true, "DELETE; -all; PRINT; -warnings 0;");
this->FindComponents();

13
src/SimParams.cpp
View File

@ -26,6 +26,7 @@
#include <iostream>
#include <string>
#include <string_view>
#include <vector>
using namespace poet;
@ -229,18 +230,18 @@ void SimParams::initVectorParams(RInside &R, int col_count) {
void SimParams::setDtDiffer(bool dt_differ) { simparams.dt_differ = dt_differ; }
t_simparams SimParams::getNumParams() { return this->simparams; }
t_simparams SimParams::getNumParams() const { return this->simparams; }
std::vector<int> SimParams::getDHTSignifVector() {
std::vector<int> SimParams::getDHTSignifVector() const {
return this->dht_signif_vector;
}
std::vector<std::string> SimParams::getDHTPropTypeVector() {
std::vector<std::string> SimParams::getDHTPropTypeVector() const {
return this->dht_prop_type_vector;
}
std::string SimParams::getDHTFile() { return this->dht_file; }
std::string_view SimParams::getDHTFile() const { return this->dht_file; }
std::string SimParams::getFilesim() { return this->filesim; }
std::string SimParams::getOutDir() { return this->out_dir; }
std::string_view SimParams::getFilesim() const { return this->filesim; }
std::string_view SimParams::getOutDir() const { return this->out_dir; }
std::list<std::string> SimParams::validateOptions(argh::parser cmdl) {
/* store all unknown parameters here */

54
util/data_evaluation/RFun_Eval.R
View File

@ -1,10 +1,20 @@
## Simple library of functions to assess and visualize the results of the coupled simulations
## Time-stamp: "Last modified 2022-12-15 11:30:55 delucia"
## Time-stamp: "Last modified 2023-01-18 16:02:58 delucia"
require(RedModRphree)
require(Rmufits) ## essentially for PlotCartCellData
require(Rcpp)
curdir <- dirname(sys.frame(1)$ofile) ##path.expand(".")
print(paste("RFun_Eval.R is in ", curdir))
sourceCpp(file = paste0(curdir, "/interpret_keys.cpp"))
# Wrapper around previous sourced Rcpp function
ConvertDHTKey <- function(value) {
rcpp_key_convert(value)
}
## function which reads all simulation results in a given directory
ReadRTSims <- function(dir) {
@ -16,16 +26,16 @@ ReadRTSims <- function(dir) {
}
## function which reads all successive DHT stored in a given directory
ReadAllDHT <- function(dir) {
ReadAllDHT <- function(dir, new_scheme = TRUE) {
files_full <- list.files(dir, pattern="iter.*dht", full.names=TRUE)
files_name <- list.files(dir, pattern="iter.*dht", full.names=FALSE)
res <- lapply(files_full, ReadDHT)
res <- lapply(files_full, ReadDHT, new_scheme = new_scheme)
names(res) <- gsub(".rds","",files_name, fixed=TRUE)
return(res)
}
## function which reads one .dht file and gives a matrix
ReadDHT <- function(file) {
ReadDHT <- function(file, new_scheme = TRUE) {
conn <- file(file, "rb") ## open for reading in binary mode
if (!isSeekable(conn))
stop("Connection not seekable")
@ -46,6 +56,15 @@ ReadDHT <- function(file) {
## close connection
close(conn)
res <- matrix(buff, nrow=nrow, ncol=ncol, byrow=TRUE)
if (new_scheme) {
nkeys <- dims[1] / 8
keys <- res[, 1:nkeys]
conv <- apply(keys, 2, ConvertDHTKey)
res[, 1:nkeys] <- conv
}
return(res)
}
@ -74,27 +93,27 @@ PlotScatter <- function(sam1, sam2, which=NULL, labs=c("NO DHT", "DHT"), pch="."
##### Some metrics for relative comparison
## Using range as norm
RranRMSE <- function (pred, obs)
RranRMSE <- function(pred, obs)
sqrt(mean((pred - obs)^2))/abs(max(pred) - min(pred))
## Using max val as norm
RmaxRMSE <- function (pred, obs)
RmaxRMSE <- function(pred, obs)
sqrt(mean((pred - obs)^2)/abs(max(pred)))
## Using sd as norm
RsdRMSE <- function (pred, obs)
RsdRMSE <- function(pred, obs)
sqrt(mean((pred - obs)^2))/sd(pred)
## Using mean as norm
RmeanRMSE <- function (pred, obs)
RmeanRMSE <- function(pred, obs)
sqrt(mean((pred - obs)^2))/mean(pred)
## Using mean as norm
RAEmax <- function (pred, obs)
RAEmax <- function(pred, obs)
mean(abs(pred - obs))/max(pred)
## Max absolute error
MAE <- function (pred, obs)
MAE <- function(pred, obs)
max(abs(pred - obs))
## workhorse function for ComputeErrors and its use with mapply
@ -161,7 +180,7 @@ ExportToParaview <- function(vtu, nameout, results) {
## "breaks" for color coding of 2D simulations
Plot2DCellData <- function (data, grid, nx, ny, contour = TRUE,
nlevels = 12, breaks, palette = "heat.colors",
rev.palette = TRUE, scale = TRUE, ...) {
rev.palette = TRUE, scale = TRUE, plot.axes=TRUE, ...) {
if (!missing(grid)) {
xc <- unique(sort(grid$cell$XCOORD))
yc <- unique(sort(grid$cell$YCOORD))
@ -190,11 +209,14 @@ Plot2DCellData <- function (data, grid, nx, ny, contour = TRUE,
1))
}
par(las = 1, mar = c(5, 5, 3, 1))
image(xc, yc, pp, xlab = "X [m]", ylab = "Y[m]", las = 1,
asp = 1, breaks = breaks, col = colors, axes = FALSE,
...)
axis(1)
axis(2)
image(xc, yc, pp, xlab = "X [m]", ylab = "Y[m]", las = 1, asp = 1,
breaks = breaks, col = colors, axes = FALSE, ann=plot.axes,
...)
if (plot.axes) {
axis(1)
axis(2)
}
if (contour)
contour(unique(sort(xc)), unique(sort(yc)), pp, breaks = breaks,
add = TRUE)

37
util/data_evaluation/interpret_keys.cpp Normal file
View File

@ -0,0 +1,37 @@
/* This file is intended to be used as input file for 'Rcpp::sourceCPP()'.
*
* The provided function will translate our key data structure back into human
* readable double values also interpretable by R or other languages.
*/
#include <Rcpp.h>
#include <cmath>
#include <cstdint>
#include <vector>
using DHT_Keyelement = struct keyelem {
std::int8_t exp : 8;
std::int64_t significant : 56;
};
using namespace Rcpp;
// [[Rcpp::export]]
std::vector<double> rcpp_key_convert(std::vector<double> input) {
std::vector<double> output;
output.reserve(input.size());
for (const double &value : input) {
DHT_Keyelement currKeyelement = *((DHT_Keyelement *)&value);
double normalize =
((std::int32_t)-std::log10(std::fabs(currKeyelement.significant))) +
currKeyelement.exp;
output.push_back(currKeyelement.significant * std::pow(10., normalize));
}
return output;
}