Compare commits

...

11 Commits

Author SHA1 Message Date
Daos Test
7c99862388 Added benchmark from Max 2023-09-06 00:08:45 +02:00
Max Luebke
1339cf80ae Update submodules 2023-08-04 13:42:00 +02:00
Daos Test
477e4211f5 "Merged" V0.x into nico-daos
It Does Compile
2023-08-03 16:42:30 +02:00
Max Lübke
3f8d6cd7a8 Revert "Merge branch 'merge-ma-changes' into nico-daos"
This reverts merge request !18
2023-07-27 15:23:05 +02:00
Max Luebke
6084c8ac80 Merge branch 'merge-ma-changes' into nico-daos 2023-07-27 15:22:22 +02:00
Max Luebke
afe1d28677 perf: handling of rounded input keys (DHT) improved
data: moved benchmarks from `data` to `bench`

build: bump PhreeqcRM to latest version

fix: removed unnecessary functions from ChemistryModule

fix: remove misses from DHT counter, as they can be deduced from count
of hits
2023-07-27 15:10:38 +02:00
Daos Test
0e69037273 Added KV-API 2023-07-09 22:01:01 +02:00
Daos Test
d54ca353bb Modified CMAKE for different DAOS implementations 2023-07-07 15:27:48 +02:00
nsauerbrei
d6b09a12d7 Bug fixes - Functional State 2023-06-21 23:36:00 +02:00
nsauerbrei
9014948805 DAOS implementation 2023-06-19 23:17:04 +02:00
nsauerbrei
e81e7ec2fa Readme Mod 2023-06-19 13:49:43 +02:00
107 changed files with 14125 additions and 1354 deletions

6
.vscode/settings.json vendored Normal file
View File

@ -0,0 +1,6 @@
{
"files.associations": {
"inttypes.h": "c",
"stdio.h": "c"
}
}

0
CMake/FindRRuntime.cmake Normal file → Executable file
View File

2
CMake/POET_Scripts.cmake Normal file → Executable file
View File

@ -9,7 +9,7 @@ macro(get_POET_version)
OUTPUT_VARIABLE POET_GIT_BRANCH
OUTPUT_STRIP_TRAILING_WHITESPACE)
execute_process(
COMMAND ${GIT_EXECUTABLE} describe --always
COMMAND ${GIT_EXECUTABLE} describe --tags --always
WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}
OUTPUT_VARIABLE POET_GIT_VERSION
OUTPUT_STRIP_TRAILING_WHITESPACE)

1
CMakeLists.txt Normal file → Executable file
View File

@ -24,7 +24,6 @@ find_package(RRuntime REQUIRED)
add_subdirectory(src)
add_subdirectory(R_lib)
add_subdirectory(data)
add_subdirectory(app)
add_subdirectory(bench)

0
LICENSE Normal file → Executable file
View File

201
README.md
View File

@ -1,200 +1,43 @@
<!--
Time-stamp: "Last modified 2023-07-21 12:34:43 mluebke"
Time-stamp: "Last modified 2023-01-19 12:06:10 delucia"
-->
# POET
# POET Daos
[POET](https://doi.org/10.5281/zenodo.4757913) is a coupled reactive transport
simulator implementing a parallel architecture and a fast, original MPI-based
Distributed Hash Table.
## Branch Information
![POET's Coupling Scheme](./docs/20230720_Scheme_POET_en.svg)
## Parsed code documentiation
A parsed version of POET's documentiation can be found at [Gitlab
pages](https://naaice.git-pages.gfz-potsdam.de/poet).
## External Libraries
The following external header library is shipped with POET:
- **argh** - https://github.com/adishavit/argh (BSD license)
- **PhreeqcRM** with patches from GFZ -
https://www.usgs.gov/software/phreeqc-version-3 -
https://git.gfz-potsdam.de/mluebke/phreeqcrm-gfz
- **tug** - https://git.gfz-potsdam.de/sec34/tug
This branch replaces POET's MPI-based key-value store with the one provided by
the DAOS-API, which is meant to be run in an environment which is connected with
a DAOS server.
## Installation
### Requirements
To install this version of POET please follow the instructions in the main branche.
The DAOS library must be build on the system.
To compile POET you need several software to be installed:
- C/C++ compiler (tested with GCC)
- MPI-Implementation (tested with OpenMPI and MVAPICH)
- R language and environment
- CMake 3.9+
- *optional*: `doxygen` with `dot` bindings for documentiation
The following R libraries must then be installed, which will get the
needed dependencies automatically:
- [Rcpp](https://cran.r-project.org/web/packages/Rcpp/index.html)
- [RInside](https://cran.r-project.org/web/packages/RInside/index.html)
### Compiling source code
The generation of makefiles is done with CMake. You should be able to generate
Makefiles by running:
To build POET with DAOS as the key-value store either the POET_DAOSKV or POET_DAOSOB flag needs to be true,
where POET_DAOSKV uses DAOS KV API, and POET_DAOSOB uses DAOS object API. For example:
```sh
mkdir build && cd build
cmake ..
cmake -DPOET_DAOSOB=True -DCMAKE_INSTALL_PREFIX=/home/dtest/poet ..
```
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 everything went well you'll find the executable at
`build/app/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:
- **POET_DHT_Debug**=_boolean_ - toggles the output of detailed statistics about
DHT usage. Defaults to _OFF_.
- **POET_ENABLE_TESTING**=_boolean_ - enables small set of unit tests (more to
come). Defaults to _OFF_.
### 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:
```sh
# start R environment
$ R
# install R dependencies
> install.packages(c("Rcpp", "RInside"))
> q(save="no")
# cd into POET project root
$ cd <POET_dir>
# Build process
$ mkdir build && cd build
$ cmake -DCMAKE_INSTALL_PREFIX=/home/<user>/poet ..
$ 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.
The correspondending directory tree would look like this:
```sh
poet
├── bin
│ └── poet
├── R_lib
│ └── kin_r_library.R
└── share
└── poet
├── bench
│ ├── dolo_diffu_inner_large.R
│ ├── dolo_diffu_inner.R
│ └── dolo_inner.pqi
└── examples
├── dol.pqi
├── phreeqc_kin.dat
├── SimDol1D_diffu.R
└── SimDol2D_diffu.R
```
The R libraries will be loaded at runtime and the paths are hardcoded
absolute paths inside `poet.cpp`. So, if you consider to move
`bin/poet` either change paths of the R source files and recompile
POET or also move `R_lib/*` relative to the binary.
## Running
Run POET by `mpirun ./poet <OPTIONS> <SIMFILE> <OUTPUT_DIRECTORY>`
where:
- **OPTIONS** - runtime parameters (explained below)
- **SIMFILE** - simulation described as R script (e.g.
`<POET_INSTALL_DIR>/share/examples/SimDol2D_diffu.R`)
- **OUTPUT_DIRECTORY** - path, where all output of POET should be stored
### Runtime options
The following parameters can be set:
| Option | Value | Description |
|--------------------------|--------------|--------------------------------------------------------------------------------------------------------------------------|
| **--work-package-size=** | _1..n_ | size of work packages (defaults to _5_) |
| **--ignore-result** | | disables store of simulation resuls |
| **--dht** | | enabling DHT usage (defaults to _OFF_) |
| **--dht-signif=** | _1..n_ | set rounding to number of significant digits (defaults to _5_) (it is recommended to use `signif_vec` in R input script) |
| **--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-snaps=** | _0-2_ | disable or enable storage of DHT snapshots |
| **--dht-file=** | `<SNAPSHOT>` | initializes DHT with the given snapshot file |
#### 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`.
#### Additions to `dht-snaps`
Following values can be set:
- _0_ = snapshots are disabled
- _1_ = only stores snapshot at the end of the simulation with name
`<OUTPUT_DIRECTORY>.dht`
- _2_ = stores snapshot at the end and after each iteration iteration
snapshot files are stored in `<DIRECTORY>/iter<n>.dht`
### Example: Running from scratch
We will continue the above example and start a simulation with
`SimDol2D_diffu.R`. As transport a simple fixed-coefficient diffusion is used.
It's a 2D, 100x100 grid, simulating 10 time steps. To start the simulation with
4 processes `cd` into your previously installed POET-dir
`<POET_INSTALL_DIR>/bin` and run:
To run POET with DAOS the --dht (WIP) parameter needs to be set, furthermore,
-x DAOS_POOL needs to be given as well, with DAOS_POOL being an environment variable with a running pool name,
for example:
```sh
mpirun -n 4 ./poet ../share/poet/examples/SimDol2D_diffu.R output
export DAOS_POOL="<pool_name>"
mpirun -n 4 -x DAOS_POOL ./poet --dht ../bench/dolo_diffu_inner/dolo_diffu_inner.R output
```
After a finished simulation all data generated by POET will be found
in the directory `output`.
## Differences to the main branch
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. Also, after each iteration a DHT snapshot shall be
produced. This is done by appending the `--dht-snaps=<value>` option. The
resulting call would look like this:
- src/ChemistryModule/CMakeLists.txt was modified to compile DAOS code
- src/ChemistryModule/DaosKeyValue.c implements the underlying framework to connect and disconnect from the DAOS server, as well as call the read and write operations
- src/ChemistryModule/DHT_Wrapper.cpp was slightly modified to use DaosKeyValue.c instead of DHT.c
```sh
mpirun -n 4 ./poet --dht --dht-snaps=2 ../share/poet/examples/SimDol2D_diffu.R output
```
## About the usage of MPI_Wtime()
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).
- include/poet/ was with the header file of DaosKeyValue extend and DHT_Wrapper.hpp was modified as well

0
R_lib/CMakeLists.txt Normal file → Executable file
View File

47
R_lib/kin_r_library.R Normal file → Executable file
View File

@ -64,28 +64,33 @@ master_init <- function(setup) {
## calculated time step if store_result is TRUE and increments the
## iteration counter
master_iteration_end <- function(setup) {
iter <- setup$iter
## MDL Write on disk state_T and state_C after every iteration
## comprised in setup$out_save
if (setup$store_result) {
if (iter %in% setup$out_save) {
nameout <- paste0(fileout, "/iter_", sprintf("%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, ">")
iter <- setup$iter
## max digits for iterations
dgts <- as.integer(ceiling(log10(setup$iterations + 1)))
## string format to use in sprintf
fmt <- paste0("%0", dgts, "d")
## Write on disk state_T and state_C after every iteration
## comprised in setup$out_save
if (setup$store_result) {
if (iter %in% setup$out_save) {
nameout <- paste0(fileout, "/iter_", sprintf(fmt=fmt, iter), ".rds")
info <- list(
tr_req_dt = as.integer(setup$req_dt)
## tr_allow_dt = setup$allowed_dt,
## tr_inniter = as.integer(setup$inniter)
)
saveRDS(list(
T = setup$state_T, C = setup$state_C,
simtime = as.integer(setup$simtime),
tr_info = info
), file = nameout)
msgm("results stored in <", nameout, ">")
}
}
}
msgm("done iteration", iter, "/", setup$maxiter)
setup$iter <- setup$iter + 1
return(setup)
msgm("done iteration", iter, "/", setup$maxiter)
setup$iter <- setup$iter + 1
return(setup)
}
## function for the workers to compute chemistry through PHREEQC

0
R_lib/kin_r_library_warnings.R Normal file → Executable file
View File

0
app/CMakeLists.txt Normal file → Executable file
View File

87
app/poet.cpp Normal file → Executable file
View File

@ -70,7 +70,7 @@ void writeFieldsToR(RInside &R, const Field &trans, const Field &chem) {
void set_chem_parameters(poet::ChemistryModule &chem, uint32_t wp_size,
const std::string &database_path) {
chem.SetErrorHandlerMode(1);
chem.SetComponentH2O(true);
chem.SetComponentH2O(false);
chem.SetRebalanceFraction(0.5);
chem.SetRebalanceByCell(true);
chem.UseSolutionDensityVolume(false);
@ -112,7 +112,6 @@ void set_chem_parameters(poet::ChemistryModule &chem, uint32_t wp_size,
}
inline double RunMasterLoop(SimParams &params, RInside &R,
ChemistryParams &chem_params,
const GridParams &g_params, uint32_t nxyz_master) {
DiffusionParams d_params{R};
@ -123,10 +122,11 @@ inline double RunMasterLoop(SimParams &params, RInside &R,
double sim_time = .0;
ChemistryModule chem(nxyz_master, params.getNumParams().wp_size,
MPI_COMM_WORLD);
set_chem_parameters(chem, nxyz_master, chem_params.database_path);
chem.RunInitFile(chem_params.input_script);
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());
@ -135,24 +135,11 @@ inline double RunMasterLoop(SimParams &params, RInside &R,
chem.setProgressBarPrintout(true);
}
if (params.getNumParams().dht_enabled) {
chem.SetDHTEnabled(true, params.getNumParams().dht_size_per_process,
chem_params.dht_species);
if (!chem_params.dht_signif.empty()) {
chem.SetDHTSignifVector(chem_params.dht_signif);
}
if (!params.getDHTFile().empty()) {
chem.ReadDHTFile(params.getDHTFile());
}
if (params.getNumParams().dht_snaps > 0) {
chem.SetDHTSnaps(params.getNumParams().dht_snaps, params.getOutDir());
}
}
/* SIMULATION LOOP */
double dStartTime = MPI_Wtime();
double dSimTime{0};
for (uint32_t iter = 1; iter < maxiter + 1; iter++) {
double start_t = MPI_Wtime();
uint32_t tick = 0;
// cout << "CPP: Evaluating next time step" << endl;
// R.parseEvalQ("mysetup <- master_iteration_setup(mysetup)");
@ -200,7 +187,8 @@ inline double RunMasterLoop(SimParams &params, RInside &R,
<< endl;
// MPI_Barrier(MPI_COMM_WORLD);
double end_t = MPI_Wtime();
dSimTime += end_t - start_t;
} // END SIMULATION LOOP
R.parseEvalQ("profiling <- list()");
@ -229,23 +217,36 @@ inline double RunMasterLoop(SimParams &params, RInside &R,
// R["phreeqc_count"] = phreeqc_counts;
// R.parseEvalQ("profiling$phreeqc_count <- phreeqc_count");
if (params.getNumParams().dht_enabled) {
if (params.getChemParams().use_dht) {
R["dht_hits"] = Rcpp::wrap(chem.GetWorkerDHTHits());
R.parseEvalQ("profiling$dht_hits <- dht_hits");
R["dht_miss"] = Rcpp::wrap(chem.GetWorkerDHTMiss());
R.parseEvalQ("profiling$dht_miss <- dht_miss");
R["dht_evictions"] = Rcpp::wrap(chem.GetWorkerDHTEvictions());
R.parseEvalQ("profiling$dht_evictions <- dht_evictions");
R["dht_get_time"] = Rcpp::wrap(chem.GetWorkerDHTGetTimings());
R.parseEvalQ("profiling$dht_get_time <- dht_get_time");
R["dht_fill_time"] = Rcpp::wrap(chem.GetWorkerDHTFillTimings());
R.parseEvalQ("profiling$dht_fill_time <- dht_fill_time");
if (params.getChemParams().use_interp) {
R["interp_w"] = Rcpp::wrap(chem.GetWorkerInterpolationWriteTimings());
R.parseEvalQ("profiling$interp_write <- interp_w");
R["interp_r"] = Rcpp::wrap(chem.GetWorkerInterpolationReadTimings());
R.parseEvalQ("profiling$interp_read <- interp_r");
R["interp_g"] = Rcpp::wrap(chem.GetWorkerInterpolationGatherTimings());
R.parseEvalQ("profiling$interp_gather <- interp_g");
R["interp_fc"] =
Rcpp::wrap(chem.GetWorkerInterpolationFunctionCallTimings());
R.parseEvalQ("profiling$interp_function_calls <- interp_fc");
R["interp_calls"] = Rcpp::wrap(chem.GetWorkerInterpolationCalls());
R.parseEvalQ("profiling$interp_calls <- interp_calls");
R["interp_cached"] = Rcpp::wrap(chem.GetWorkerPHTCacheHits());
R.parseEvalQ("profiling$interp_cached <- interp_cached");
}
}
chem.MasterLoopBreak();
diffusion.end();
return MPI_Wtime() - dStartTime;
return dSimTime;
}
int main(int argc, char *argv[]) {
@ -265,19 +266,23 @@ int main(int argc, char *argv[]) {
if (world_rank > 0) {
{
uint32_t c_size;
MPI_Bcast(&c_size, 1, MPI_UINT32_T, 0, MPI_COMM_WORLD);
char *buffer = new char[c_size + 1];
MPI_Bcast(buffer, c_size, MPI_CHAR, 0, MPI_COMM_WORLD);
buffer[c_size] = '\0';
SimParams params(world_rank, world_size);
int pret = params.parseFromCmdl(argv, RInsidePOET::getInstance());
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);
set_chem_parameters(worker, worker.GetWPSize(), std::string(buffer));
delete[] buffer;
ChemistryModule worker = poet::ChemistryModule::createWorker(
MPI_COMM_WORLD, params.getChemParams());
set_chem_parameters(worker, worker.GetWPSize(),
params.getChemParams().database_path);
worker.WorkerLoop();
}
@ -331,17 +336,9 @@ int main(int argc, char *argv[]) {
// MPI_Barrier(MPI_COMM_WORLD);
poet::ChemistryParams chem_params(R);
/* THIS IS EXECUTED BY THE MASTER */
std::string db_path = chem_params.database_path;
uint32_t c_size = db_path.size();
MPI_Bcast(&c_size, 1, MPI_UINT32_T, 0, MPI_COMM_WORLD);
MPI_Bcast(db_path.data(), c_size, MPI_CHAR, 0, MPI_COMM_WORLD);
uint32_t nxyz_master = (world_size == 1 ? g_params.total_n : 1);
dSimTime = RunMasterLoop(params, R, chem_params, g_params, nxyz_master);
dSimTime = RunMasterLoop(params, R, g_params, nxyz_master);
cout << "CPP: finished simulation loop" << endl;

0
app/poet.h.in Normal file → Executable file
View File

3
bench/CMakeLists.txt Normal file → Executable file
View File

@ -1,3 +1,4 @@
add_subdirectory(dolo_diffu_inner)
add_subdirectory(dolo)
add_subdirectory(surfex)
add_subdirectory(barite)
add_subdirectory(DAOSBenchmark)

View File

@ -0,0 +1,7 @@
install(FILES
dolo_diffu_edge.R
dolo_inner.pqi
phreeqc_kin.dat
DESTINATION
share/poet/bench/daos
)

View File

@ -0,0 +1,205 @@
## Time-stamp: "Last modified 2023-09-05 14:42:20 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, 1, 1)
# 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) ##
#################################################################
## # 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(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
# pht_species = pht_species
)
#################################################################
## Section 4 ##
## Putting all those things together ##
#################################################################
iterations <- 1500
dt <- 500
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))
)

View File

@ -2,7 +2,7 @@
### SURFACE and with the RATES for Calcite and Dolomite to use with
### RedModRphree
### Time-stamp: "Last modified 2018-05-06 14:36:23 delucia"
### Time-stamp: "Last modified 2023-05-23 10:35:56 mluebke"
# PHREEQC.DAT for calculating pressure dependence of reactions, with
# molal volumina of aqueous species and of minerals, and
@ -1276,9 +1276,9 @@ Calcite
110 logK25=-5.81
120 mech_c=(10^logK25)*(e^(-Ea/R*deltaT))
130 rate=mech_a+mech_c
## 140 IF SI("Calcite")<0 then moles=parm(1)*rate*(1-SR("Calcite")) # dissolution
140 IF (SI("Calcite")<0 AND M>0) then moles=parm(1)*rate*(1-SR("Calcite")) # dissolution
## 145 IF SI("Calcite")>0 then moles=parm(1)*M*rate*(-1+SR("Calcite")) # precipitation
150 moles=parm(1)*rate*(1-SR("Calcite")) # precipitation
## 150 moles=parm(1)*rate*(1-SR("Calcite")) # precipitation
200 save moles*time
-end

1
bench/barite/CMakeLists.txt Normal file → Executable file
View File

@ -1,5 +1,6 @@
install(FILES
barite.R
barite_interp_eval.R
barite.pqi
db_barite.dat
DESTINATION

24
bench/barite/barite.R Normal file → Executable file
View File

@ -1,4 +1,4 @@
## Time-stamp: "Last modified 2023-04-24 16:51:23 mluebke"
## Time-stamp: "Last modified 2023-08-02 13:59:22 mluebke"
database <- normalizePath("../share/poet/bench/barite/db_barite.dat")
input_script <- normalizePath("../share/poet/bench/barite/barite.pqi")
@ -28,11 +28,7 @@ init_cell <- list(
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
type = types[1]
)
@ -115,16 +111,20 @@ diffusion <- list(
## # Needed when using DHT
dht_species <- names(init_diffu)
#dht_signif <- rep(15, length(dht_species))
dht_signif <- c(10, 10, 3, 5, 5, 5, 5)
dht_species <- c(
"H" = 10,
"O" = 10,
"Charge" = 3,
"Ba" = 5,
"Cl" = 5,
"S(6)" = 5,
"Sr" = 5
)
chemistry <- list(
database = database,
input_script = input_script,
dht_species = dht_species,
dht_signif = dht_signif
dht_species = dht_species
)
#################################################################

0
bench/barite/barite.pqi Normal file → Executable file
View File

151
bench/barite/barite_interp_eval.R Executable file
View File

@ -0,0 +1,151 @@
## Time-stamp: "Last modified 2023-07-21 15:04:49 mluebke"
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.217E-09,
"Ba" = 1.E-10,
"Cl" = 2.E-10,
"S" = 6.205E-4,
"Sr" = 6.205E-4,
"Barite" = 0.001,
"Celestite" = 1
)
grid <- list(
n_cells = c(n, m),
s_cells = c(n / 10, m / 10),
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,
"H" = 0.00000028904,
# "O" = 55.5087,
"O" = 0.000000165205,
# "Charge" = -1.217E-09,
"Charge" = -3.337E-08,
"Ba" = 1.E-10,
"Cl" = 1.E-10,
"S(6)" = 6.205E-4,
"Sr" = 6.205E-4
)
injection_diff <- list(
list(
# "H" = 111.0124,
"H" = 0.0000002890408,
# "O" = 55.50622,
"O" = 0.00002014464,
# "Charge" = -3.337E-08,
"Charge" = -3.337000004885E-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, floor(n / 2), floor(m / 2))
## l2 = c(2,80,80),
## l3 = c(2,60,80)
)
boundary <- list(
# "N" = rep(1, n),
"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, 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" = 10,
"O" = 10,
"Charge" = 3,
"Ba" = 5,
"Cl" = 5,
"S(6)" = 5,
"Sr" = 5
)
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)
)

0
bench/barite/db_barite.dat Normal file → Executable file
View File

View File

@ -2,6 +2,8 @@ install(FILES
dolo_diffu_inner.R
dolo_diffu_inner_large.R
dolo_inner.pqi
dolo_interp_long.R
phreeqc_kin.dat
DESTINATION
share/poet/bench/dolo
)

0
bench/dolo_diffu_inner/Eval.R → bench/dolo/Eval.R Normal file → Executable file
View File

0
data/dol.pqi → bench/dolo/dol.pqi Normal file → Executable file
View File

View File

@ -1,6 +1,6 @@
## Time-stamp: "Last modified 2023-04-24 15:12:02 mluebke"
## Time-stamp: "Last modified 2023-08-02 13:59:02 mluebke"
database <- normalizePath("../share/poet/examples/phreeqc_kin.dat")
database <- normalizePath("../share/poet/bench/dolo/phreeqc_kin.dat")
input_script <- normalizePath("../share/poet/bench/dolo/dolo_inner.pqi")
#################################################################
@ -29,11 +29,7 @@ init_cell <- list(
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
type = types[1]
)
@ -44,8 +40,8 @@ grid <- list(
## initial conditions
init_diffu <- list(
"H" = 0.000211313883539788,
"O" = 0.00398302904424952,
"H" = 110.683,
"O" = 55.3413,
"Charge" = -5.0822e-19,
"C(4)" = 1.2279E-4,
"Ca" = 1.2279E-4,
@ -67,8 +63,8 @@ alpha_diffu <- c(
## list of boundary conditions/inner nodes
vecinj_diffu <- list(
list(
"H" = 0.0001540445,
"O" = 0.002148006,
"H" = 110.683,
"O" = 55.3413,
"Charge" = 1.90431e-16,
"C(4)" = 0,
"Ca" = 0,
@ -76,8 +72,8 @@ vecinj_diffu <- list(
"Mg" = 0.001
),
list(
"H" = 0.0001610193,
"O" = 0.002386934,
"H" = 110.683,
"O" = 55.3413,
"Charge" = 1.90431e-16,
"C(4)" = 0,
"Ca" = 0.0,
@ -120,16 +116,21 @@ diffusion <- list(
## # Needed when using DHT
dht_species <- c("H", "O", "Charge", "C(4)", "Ca", "Cl", "Mg", "Calcite",
"Dolomite")
#dht_signif <- rep(15, length(dht_species))
dht_signif <- c(10, 10, 3, 5, 5, 5, 5, 5, 5)
dht_species <- c("H" = 10,
"O" = 10,
"Charge" = 3,
"C(4)" = 5,
"Ca" = 5,
"Cl" = 5,
"Mg" = 5,
"Calcite" = 5,
"Dolomite" =5
)
chemistry <- list(
database = database,
input_script = input_script,
dht_species = dht_species,
dht_signif = dht_signif
dht_species = dht_species
)
#################################################################

View File

@ -1,6 +1,6 @@
## Time-stamp: "Last modified 2023-04-24 15:43:26 mluebke"
## Time-stamp: "Last modified 2023-08-02 13:59:12 mluebke"
database <- normalizePath("../share/poet/examples/phreeqc_kin.dat")
database <- normalizePath("../share/poet/bench/dolo/phreeqc_kin.dat")
input_script <- normalizePath("../share/poet/bench/dolo/dolo_inner.pqi")
#################################################################
@ -29,11 +29,7 @@ init_cell <- list(
grid <- list(
n_cells = c(n, m),
s_cells = c(2, 1),
type = types[1],
init_cell = as.data.frame(init_cell, check.names = FALSE),
props = names(init_cell),
database = database,
input_script = input_script
type = types[1]
)
@ -118,18 +114,22 @@ diffusion <- list(
## Chemistry module (Phreeqc) ##
#################################################################
## # Needed when using DHT
dht_species <- c("H", "O", "Charge", "C(4)", "Ca", "Cl", "Mg", "Calcite",
"Dolomite")
#dht_signif <- rep(15, length(dht_species))
dht_signif <- c(10, 10, 3, 5, 5, 5, 5, 5, 5)
dht_species <- c("H" = 10,
"O" = 10,
"Charge" = 3,
"C(4)" = 5,
"Ca" = 5,
"Cl" = 5,
"Mg" = 5,
"Calcite" = 5,
"Dolomite" =5
)
chemistry <- list(
database = database,
input_script = input_script,
dht_species = dht_species,
dht_signif = dht_signif
dht_species = dht_species
)
#################################################################

28
bench/dolo/dolo_inner.pqi Executable file
View File

@ -0,0 +1,28 @@
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

167
bench/dolo/dolo_interp_long.R Executable file
View File

@ -0,0 +1,167 @@
## Time-stamp: "Last modified 2023-08-02 13:47:06 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
#)
chemistry <- list(
database = database,
input_script = input_script,
dht_species = dht_species
# 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))
)

1307
bench/dolo/phreeqc_kin.dat Executable file

File diff suppressed because it is too large Load Diff

0
bench/surfex/CMakeLists.txt Normal file → Executable file
View File

0
bench/surfex/ExBase.pqi Normal file → Executable file
View File

0
bench/surfex/SMILE_2021_11_01_TH.dat Normal file → Executable file
View File

0
bench/surfex/SurfExBase.pqi Normal file → Executable file
View File

8
bench/surfex/ex.R Normal file → Executable file
View File

@ -1,4 +1,4 @@
## Time-stamp: "Last modified 2023-04-17 12:29:27 mluebke"
## Time-stamp: "Last modified 2023-08-02 13:59:35 mluebke"
database <- normalizePath("./SMILE_2021_11_01_TH.dat")
input_script <- normalizePath("./ExBase.pqi")
@ -40,11 +40,7 @@ init_cell <- list(H = 1.476571028625e-01,
grid <- list(
n_cells = c(n, m),
s_cells = c(1, 1),
type = "scratch",
init_cell = as.data.frame(init_cell, check.names = FALSE),
props = names(init_cell),
database = database,
input_script = input_script
type = "scratch"
)

8
bench/surfex/surfex.R Normal file → Executable file
View File

@ -1,4 +1,4 @@
## Time-stamp: "Last modified 2023-04-17 15:48:21 mluebke"
## Time-stamp: "Last modified 2023-08-02 13:59:44 mluebke"
database <- normalizePath("../share/poet/bench/surfex/SMILE_2021_11_01_TH.dat")
input_script <- normalizePath("../share/poet/bench/surfex/SurfExBase.pqi")
@ -40,11 +40,7 @@ init_cell <- list(H = 1.476571028625e-01,
grid <- list(
n_cells = c(n, m),
s_cells = c(1, 1),
type = "scratch",
init_cell = as.data.frame(init_cell, check.names = FALSE),
props = names(init_cell),
database = database,
input_script = input_script
type = "scratch"
)

BIN
bin/poet Executable file

Binary file not shown.

View File

@ -1,8 +0,0 @@
install(
FILES
SimDol2D_diffu.R
SimDol1D_diffu.R
phreeqc_kin.dat
dol.pqi
DESTINATION
share/poet/examples)

View File

@ -1,180 +0,0 @@
#################################################################
## Section 1 ##
## Grid initialization ##
#################################################################
n <- 5
m <- 5
types <- c("scratch", "phreeqc", "rds")
# initsim <- c("SELECTED_OUTPUT", "-high_precision true",
# "-reset false",
# "USER_PUNCH",
# "-head C Ca Cl Mg pH pe O2g Calcite Dolomite",
# "10 PUNCH TOT(\"C\"), TOT(\"Ca\"), TOT(\"Cl\"), TOT(\"Mg\"), -LA(\"H+\"), -LA(\"e-\"), EQUI(\"O2g\"), EQUI(\"Calcite\"), EQUI(\"Dolomite\")",
# "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.6788 10.0",
# "Calcite 0.0 2.07E-3",
# "Dolomite 0.0 0.0",
# "END")
# needed if init type is set to "scratch"
# prop <- c("C", "Ca", "Cl", "Mg", "pH", "pe", "O2g", "Calcite", "Dolomite")
init_cell <- list(
"C" = 1.2279E-4,
"Ca" = 1.2279E-4,
"Cl" = 0,
"Mg" = 0,
"pH" = 9.91,
"pe" = 4,
"O2g" = 10,
"Calcite" = 2.07e-4,
"Dolomite" = 0
)
grid <- list(
n_cells = n,
s_cells = n,
type = types[1],
init_cell = as.data.frame(init_cell),
props = names(init_cell),
init_script = NULL
)
##################################################################
## Section 2 ##
## Diffusion parameters and boundary conditions ##
##################################################################
init_diffu <- c(
"C" = 1.2279E-4,
"Ca" = 1.2279E-4,
"Cl" = 0,
"Mg" = 0,
"pe" = 4,
"pH" = 7
)
alpha_diffu <- c(
"C" = 1E-4,
"Ca" = 1E-4,
"Cl" = 1E-4,
"Mg" = 1E-4,
"pe" = 1E-4,
"pH" = 1E-4
)
vecinj_diffu <- list(
list(
"C" = 0,
"Ca" = 0,
"Cl" = 0.002,
"Mg" = 0.001,
"pe" = 4,
"pH" = 9.91
)
)
boundary <- list(
"E" = 0,
"W" = 1
)
diffu_list <- names(alpha_diffu)
diffusion <- list(
init = init_diffu,
vecinj = do.call(rbind.data.frame, vecinj_diffu),
vecinj_index = boundary,
alpha = alpha_diffu
)
##################################################################
## Section 3 ##
## Phreeqc simulation ##
##################################################################
db <- RPhreeFile(system.file("extdata", "phreeqc_kin.dat",
package = "RedModRphree"
), is.db = TRUE)
phreeqc::phrLoadDatabaseString(db)
# NOTE: This won't be needed in the future either. Could also be done in a. pqi
# file
base <- c(
"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",
"Mg 0.001",
"Cl 0.002",
"PURE 1",
"O2g -0.1675 10",
"KINETICS 1",
"-steps 100",
"-step_divide 100",
"-bad_step_max 2000",
"Calcite", "-m 0.000207",
"-parms 0.0032",
"Dolomite",
"-m 0.0",
"-parms 0.00032",
"END"
)
selout <- c(
"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"
)
# Needed when using DHT
signif_vector <- c(7, 7, 7, 7, 7, 7, 7, 5, 5)
prop_type <- c("act", "act", "act", "act", "logact", "logact", "ignore", "act", "act")
prop <- names(init_cell)
iterations <- 500
setup <- list(
# bound = myboundmat,
base = base,
first = selout,
# initsim = initsim,
# Cf = 1,
grid = grid,
diffusion = diffusion,
prop = prop,
immobile = c(7, 8, 9),
kin = c(8, 9),
ann = list(O2g = -0.1675),
# phase = "FLUX1",
# density = "DEN1",
reduce = FALSE,
# snapshots = demodir, ## directory where we will read MUFITS SUM files
# gridfile = paste0(demodir, "/d2ascii.run.GRID.SUM")
# init = init,
# vecinj = vecinj,
# cinj = c(0,1),
# boundary = boundary,
# injections = FALSE,
iterations = iterations,
timesteps = rep(10, iterations)
)

View File

@ -1,213 +0,0 @@
database <- normalizePath("../share/poet/examples/phreeqc_kin.dat")
input_script <- normalizePath("../share/poet/examples/dol.pqi")
#################################################################
## Section 1 ##
## Grid initialization ##
#################################################################
n <- 100
m <- 100
types <- c("scratch", "phreeqc", "rds")
# initsim <- c("SELECTED_OUTPUT", "-high_precision true",
# "-reset false",
# "USER_PUNCH",
# "-head C Ca Cl Mg pH pe O2g Calcite Dolomite",
# "10 PUNCH TOT(\"C\"), TOT(\"Ca\"), TOT(\"Cl\"), TOT(\"Mg\"), -LA(\"H+\"), -LA(\"e-\"), EQUI(\"O2g\"), EQUI(\"Calcite\"), EQUI(\"Dolomite\")",
# "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.6788 10.0",
# "Calcite 0.0 2.07E-3",
# "Dolomite 0.0 0.0",
# "END")
# needed if init type is set to "scratch"
# prop <- c("C", "Ca", "Cl", "Mg", "pH", "pe", "O2g", "Calcite", "Dolomite")
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(n, m),
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 ##
##################################################################
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
)
alpha_diffu <- c(
"H" = 1E-4,
"O" = 1E-4,
"Charge" = 1E-4,
"C" = 1E-4,
"Ca" = 1E-4,
"Cl" = 1E-4,
"Mg" = 1E-4
)
vecinj_diffu <- list(
list(
"H" = 110.683,
"O" = 55.3413,
"Charge" = 1.90431e-16,
"C" = 0,
"Ca" = 0,
"Cl" = 0.002,
"Mg" = 0.001
)
)
#inner_index <- c(5, 15, 25)
#inner_vecinj_index <- rep(1, 3)
#
#vecinj_inner <- cbind(inner_index, inner_vecinj_index)
vecinj_inner <- list(
l1 = c(1,2,2)
)
boundary <- list(
"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 ##
## Chemitry module (Phreeqc) ##
#################################################################
# db <- RPhreeFile(system.file("extdata", "phreeqc_kin.dat",
# package = "RedModRphree"
# ), is.db = TRUE)
#
# phreeqc::phrLoadDatabaseString(db)
# NOTE: This won't be needed in the future either. Could also be done in a. pqi
# file
base <- c(
"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",
"Mg 0.001",
"Cl 0.002",
"PURE 1",
"O2g -0.1675 10",
"KINETICS 1",
"-steps 100",
"-step_divide 100",
"-bad_step_max 2000",
"Calcite", "-m 0.000207",
"-parms 0.0032",
"Dolomite",
"-m 0.0",
"-parms 0.00032",
"END"
)
selout <- c(
"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"
)
# Needed when using DHT
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(
database = database,
input_script = input_script
)
#################################################################
## Section 4 ##
## Putting all those things together ##
#################################################################
iterations <- 10
setup <- list(
# bound = myboundmat,
base = base,
first = selout,
# initsim = initsim,
# Cf = 1,
grid = grid,
diffusion = diffusion,
chemistry = chemistry,
prop = prop,
immobile = c(7, 8, 9),
kin = c(8, 9),
ann = list(O2g = -0.1675),
# phase = "FLUX1",
# density = "DEN1",
reduce = FALSE,
# snapshots = demodir, ## directory where we will read MUFITS SUM files
# gridfile = paste0(demodir, "/d2ascii.run.GRID.SUM")
# init = init,
# vecinj = vecinj,
# cinj = c(0,1),
# boundary = boundary,
# injections = FALSE,
iterations = iterations,
timesteps = rep(10, iterations)
)

View File

@ -1,13 +0,0 @@
phreeqc:
RebalanceByCell: True
RebalanceFraction: 0.5
SolutionDensityVolume: False
PartitionUZSolids: False
Units:
Solution: "mg/L"
PPassamblage: "mol/L"
Exchange: "mol/L"
Surface: "mol/L"
GasPhase: "mol/L"
SSassemblage: "mol/L"
Kinetics: "mol/L"

0
docs/20221216_Scheme_PORT_en.svg Normal file → Executable file
View File

Before

Width:  |  Height:  |  Size: 73 KiB

After

Width:  |  Height:  |  Size: 73 KiB

0
docs/20230720_Scheme_POET_en.svg Normal file → Executable file
View File

Before

Width:  |  Height:  |  Size: 73 KiB

After

Width:  |  Height:  |  Size: 73 KiB

0
docs/CMakeLists.txt Normal file → Executable file
View File

21
docs/Input_Scripts.md Normal file → Executable file
View File

@ -13,10 +13,6 @@ by POET.
| `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* |
| `init_cell` | Data Frame | Containing all exactly one value per species to initialize the field. |
| `props` | String Vector | Names of all species |
| `database` | String | Path to Phreeqc database |
| `input_script` | String | Path to Phreeqc initial script |
## Diffusion parameters
@ -70,10 +66,12 @@ vecinj_index <- list(
## Chemistry parameters
| name | type | description |
|----------------|--------|-----------------------------------|
| `database` | String | Path to the Phreeqc database |
| `input_script` | String | Path the the Phreeqc input script |
| 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
@ -86,10 +84,3 @@ vecinj_index <- list(
| `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 |
### DHT setup
| name | type | description |
|-----------------|----------------|---------------------------------------------------------------------------------|
| `signif_vector` | Numeric Vector | Indicates significant digits to use for DHT rounding. Order of `props` vector. |
| `prop_type` | String Vector | Set type of species for rounding, can be left blank or set to *act* or *ignore* |

0
docs/Output.md Normal file → Executable file
View File

@ -1 +1 @@
Subproject commit 8fdfd113dcb4ad1a31705ff8ddb13a21a505bad8
Subproject commit ae7a13539fb71f270b87eb2e874fbac80bc8dda2

@ -1 +1 @@
Subproject commit ba5dc40af119da015d36d2554ecd558ef9da1eb6
Subproject commit 89f713b273cd5340b2e8169523da04c2d7ad89c9

121
include/poet/ChemistryModule.hpp Normal file → Executable file
View File

@ -1,20 +1,23 @@
// Time-stamp: "Last modified 2023-07-21 17:20:10 mluebke"
// Time-stamp: "Last modified 2023-08-01 13:13:18 mluebke"
#ifndef CHEMISTRYMODULE_H_
#define CHEMISTRYMODULE_H_
#include "DHT_Wrapper.hpp"
#include "Field.hpp"
#include "Interpolation.hpp"
#include "IrmResult.h"
#include "PhreeqcRM.h"
#include "poet/DHT_Wrapper.hpp"
#include "SimParams.hpp"
#include "DataStructures.hpp"
#include <array>
#include <cstddef>
#include <cstdint>
#include <map>
#include <memory>
#include <mpi.h>
#include <string>
#include <unordered_map>
#include <vector>
namespace poet {
@ -39,7 +42,8 @@ public:
* \param wp_size Count of grid cells to fill each work package at maximum.
* \param communicator MPI communicator to distribute work in.
*/
ChemistryModule(uint32_t nxyz, uint32_t wp_size, MPI_Comm communicator);
ChemistryModule(uint32_t nxyz, uint32_t wp_size, std::uint32_t maxiter,
const ChemistryParams &chem_param, MPI_Comm communicator);
/**
* Deconstructor, which frees DHT data structure if used.
@ -86,7 +90,8 @@ public:
*
* \returns A worker instance with fixed work package size.
*/
static ChemistryModule createWorker(MPI_Comm communicator);
static ChemistryModule createWorker(MPI_Comm communicator,
const ChemistryParams &chem_param);
/**
* Default work package size.
@ -109,9 +114,9 @@ public:
* Enumerating DHT file options
*/
enum {
DHT_FILES_DISABLED = 0, //!< disabled file output
DHT_FILES_SIMEND, //!< only output of snapshot after simulation
DHT_FILES_ITEREND //!< output snapshots after each iteration
DHT_SNAPS_DISABLED = 0, //!< disabled file output
DHT_SNAPS_SIMEND, //!< only output of snapshot after simulation
DHT_SNAPS_ITEREND //!< output snapshots after each iteration
};
/**
@ -135,41 +140,6 @@ public:
*/
void MasterLoopBreak();
/**
* **Master only** Enables DHT usage for the workers.
*
* If called multiple times enabling the DHT, the current state of DHT will be
* overwritten with a new instance of the DHT.
*
* \param enable Enables or disables the usage of the DHT.
* \param size_mb Size in megabyte to allocate for the DHT if enabled.
* \param key_species Name of species to be used for key creation.
*/
void SetDHTEnabled(bool enable, uint32_t size_mb,
const std::vector<std::string> &key_species);
/**
* **Master only** Set DHT snapshots to specific mode.
*
* \param type DHT snapshot mode.
* \param out_dir Path to output DHT snapshots.
*/
void SetDHTSnaps(int type, const std::string &out_dir);
/**
* **Master only** Set the vector with significant digits to round before
* inserting into DHT.
*
* \param signif_vec Vector defining significant digit for each species. Order
* is defined by prop_type vector (ChemistryModule::GetPropNames).
*/
void SetDHTSignifVector(std::vector<uint32_t> &signif_vec);
/**
* **Master only** Load the state of the DHT from given file.
*
* \param input_file File to load data from.
*/
void ReadDHTFile(const std::string &input_file);
/**
* **Master only** Return count of grid cells.
*/
@ -231,13 +201,6 @@ public:
*/
std::vector<uint32_t> GetWorkerDHTHits() const;
/**
* **Master only** Collect and return DHT misses of all workers.
*
* \return Vector of all count of DHT misses.
*/
std::vector<uint32_t> GetWorkerDHTMiss() const;
/**
* **Master only** Collect and return DHT evictions of all workers.
*
@ -257,20 +220,43 @@ public:
*
* \param enabled True if print progressbar, false if not.
*/
void setProgressBarPrintout(bool enabled);
void setProgressBarPrintout(bool enabled) {
this->print_progessbar = enabled;
};
std::vector<uint32_t> GetWorkerInterpolationCalls() const;
std::vector<double> GetWorkerInterpolationWriteTimings() const;
std::vector<double> GetWorkerInterpolationReadTimings() const;
std::vector<double> GetWorkerInterpolationGatherTimings() const;
std::vector<double> GetWorkerInterpolationFunctionCallTimings() const;
std::vector<uint32_t> GetWorkerPHTCacheHits() const;
protected:
void initializeDHT(uint32_t size_mb,
const NamedVector<std::uint32_t> &key_species);
void setDHTSnapshots(int type, const std::string &out_dir);
void setDHTReadFile(const std::string &input_file);
void initializeInterp(std::uint32_t bucket_size, std::uint32_t size_mb,
std::uint32_t min_entries,
const NamedVector<std::uint32_t> &key_species);
enum {
CHEM_INIT,
CHEM_WP_SIZE,
CHEM_INIT_SPECIES,
CHEM_DHT_ENABLE,
CHEM_DHT_SIGNIF_VEC,
CHEM_DHT_PROP_TYPE_VEC,
CHEM_DHT_SNAPS,
CHEM_DHT_READ_FILE,
CHEM_IP_ENABLE,
CHEM_IP_MIN_ENTRIES,
CHEM_IP_SIGNIF_VEC,
CHEM_WORK_LOOP,
CHEM_PERF,
CHEM_PROGRESSBAR,
CHEM_BREAK_MAIN_LOOP
};
@ -281,11 +267,20 @@ protected:
WORKER_DHT_GET,
WORKER_DHT_FILL,
WORKER_IDLE,
WORKER_IP_WRITE,
WORKER_IP_READ,
WORKER_IP_GATHER,
WORKER_IP_FC,
WORKER_DHT_HITS,
WORKER_DHT_MISS,
WORKER_DHT_EVICTIONS
WORKER_DHT_EVICTIONS,
WORKER_PHT_CACHE_HITS,
WORKER_IP_CALLS
};
std::vector<uint32_t> interp_calls;
std::vector<uint32_t> dht_hits;
std::vector<uint32_t> dht_evictions;
struct worker_s {
double phreeqc_t = 0.;
double dht_get = 0.;
@ -340,8 +335,11 @@ protected:
void unshuffleField(const std::vector<double> &in_buffer,
uint32_t size_per_prop, uint32_t prop_count,
uint32_t wp_count, std::vector<double> &out_field);
std::vector<std::uint32_t>
parseDHTSpeciesVec(const std::vector<std::string> &species_vec) const;
std::vector<std::int32_t>
parseDHTSpeciesVec(const NamedVector<std::uint32_t> &key_species,
const std::vector<std::string> &to_compare) const;
void BCastStringVec(std::vector<std::string> &io);
int comm_size, comm_rank;
MPI_Comm group_comm;
@ -351,11 +349,14 @@ protected:
uint32_t wp_size{CHEM_DEFAULT_WP_SIZE};
bool dht_enabled{false};
int dht_snaps_type{DHT_FILES_DISABLED};
int dht_snaps_type{DHT_SNAPS_DISABLED};
std::string dht_file_out_dir;
poet::DHT_Wrapper *dht = nullptr;
bool interp_enabled{false};
std::unique_ptr<poet::InterpolationModule> interp;
static constexpr uint32_t BUFFER_OFFSET = 5;
inline void ChemBCast(void *buf, int count, MPI_Datatype datatype) const {
@ -374,7 +375,9 @@ protected:
bool print_progessbar{false};
double chem_t = 0.;
std::uint32_t file_pad;
double chem_t{0.};
uint32_t n_cells = 0;
uint32_t prop_count = 0;
@ -384,6 +387,8 @@ protected:
static constexpr int MODULE_COUNT = 5;
const ChemistryParams &params;
std::array<std::uint32_t, MODULE_COUNT> speciesPerModule{};
};
} // namespace poet

51
include/poet/DHT.h Normal file → Executable file
View File

@ -67,7 +67,7 @@
*/
typedef struct {
/** Count of writes to specific process this process did. */
int* writes_local;
int *writes_local;
/** Writes after last call of DHT_print_statistics. */
int old_writes;
/** How many read misses occur? */
@ -100,27 +100,36 @@ typedef struct {
/** Size of the MPI communicator respectively all participating processes. */
int comm_size;
/** Pointer to a hashfunction. */
uint64_t (*hash_func)(int, const void*);
uint64_t (*hash_func)(int, const void *);
/** Pre-allocated memory where a bucket can be received. */
void* recv_entry;
void *recv_entry;
/** Pre-allocated memory where a bucket to send can be stored. */
void* send_entry;
void *send_entry;
/** Allocated memory on which the MPI window was created. */
void* mem_alloc;
void *mem_alloc;
/** Count of read misses over all time. */
int read_misses;
/** Count of evictions over all time. */
int evictions;
/** Array of indeces where a bucket can be stored. */
uint64_t* index;
uint64_t *index;
/** Count of possible indeces. */
unsigned int index_count;
int (*accumulate_callback)(int, void *, int, void *);
#ifdef DHT_STATISTICS
/** Detailed statistics of the usage of the DHT. */
DHT_stats* stats;
DHT_stats *stats;
#endif
} DHT;
extern void DHT_set_accumulate_callback(DHT *table,
int (*callback_func)(int, void *, int,
void *));
extern int DHT_write_accumulate(DHT *table, const void *key, int send_size,
void *data, uint32_t *proc, uint32_t *index, int *callback_ret);
/**
* @brief Create a DHT.
*
@ -141,9 +150,9 @@ typedef struct {
* @return DHT* The returned value is the \a DHT-object which serves as a handle
* for all DHT operations. If an error occured NULL is returned.
*/
extern DHT* DHT_create(MPI_Comm comm, uint64_t size_per_process,
extern DHT *DHT_create(MPI_Comm comm, uint64_t size_per_process,
unsigned int data_size, unsigned int key_size,
uint64_t (*hash_func)(int, const void*));
uint64_t (*hash_func)(int, const void *));
/**
* @brief Write data into DHT.
@ -161,10 +170,14 @@ extern DHT* DHT_create(MPI_Comm comm, uint64_t size_per_process,
* @param table Pointer to the \a DHT-object.
* @param key Pointer to the key.
* @param data Pointer to the data.
* @param proc If not NULL, returns the process number written to.
* @param index If not NULL, returns the index of the bucket where the data was
* written to.
* @return int Returns either DHT_SUCCESS on success or correspondending error
* value on eviction or error.
*/
extern int DHT_write(DHT* table, void* key, void* data);
extern int DHT_write(DHT *table, void *key, void *data, uint32_t *proc,
uint32_t *index);
/**
* @brief Read data from DHT.
@ -187,8 +200,10 @@ extern int DHT_write(DHT* table, void* key, void* data);
* @return int Returns either DHT_SUCCESS on success or correspondending error
* value on read miss or error.
*/
extern int DHT_read(DHT* table, void* key, void* destination);
extern int DHT_read(DHT *table, const void *key, void *destination);
extern int DHT_read_location(DHT *table, uint32_t proc, uint32_t index,
void *destination);
/**
* @brief Write current state of DHT to file.
*
@ -203,7 +218,7 @@ extern int DHT_read(DHT* table, void* key, void* destination);
* @return int Returns DHT_SUCCESS on succes, DHT_FILE_IO_ERROR if file can't be
* opened/closed or DHT_WRITE_ERROR if file is not writable.
*/
extern int DHT_to_file(DHT* table, const char* filename);
extern int DHT_to_file(DHT *table, const char *filename);
/**
* @brief Read state of DHT from file.
@ -223,7 +238,7 @@ extern int DHT_to_file(DHT* table, const char* filename);
* file doesn't match expectation. This is possible if the data size or key size
* is different.
*/
extern int DHT_from_file(DHT* table, const char* filename);
extern int DHT_from_file(DHT *table, const char *filename);
/**
* @brief Free ressources of DHT.
@ -241,7 +256,7 @@ extern int DHT_from_file(DHT* table, const char* filename);
* @return int Returns either DHT_SUCCESS on success or DHT_MPI_ERROR on
* internal MPI error.
*/
extern int DHT_free(DHT* table, int* eviction_counter, int* readerror_counter);
extern int DHT_free(DHT *table, int *eviction_counter, int *readerror_counter);
/**
* @brief Prints a table with statistics about current use of DHT.
@ -267,7 +282,7 @@ extern int DHT_free(DHT* table, int* eviction_counter, int* readerror_counter);
* @return int Returns DHT_SUCCESS on success or DHT_MPI_ERROR on internal MPI
* error.
*/
extern int DHT_print_statistics(DHT* table);
extern int DHT_print_statistics(DHT *table);
/**
* @brief Determine destination rank and index.
@ -286,8 +301,8 @@ extern int DHT_print_statistics(DHT* table);
* @param index_count Count of possible indeces.
*/
static void determine_dest(uint64_t hash, int comm_size,
unsigned int table_size, unsigned int* dest_rank,
uint64_t* index, unsigned int index_count);
unsigned int table_size, unsigned int *dest_rank,
uint64_t *index, unsigned int index_count);
/**
* @brief Set the occupied flag.
@ -296,7 +311,7 @@ static void determine_dest(uint64_t hash, int comm_size,
*
* @param flag_byte First byte of a bucket.
*/
static void set_flag(char* flag_byte);
static void set_flag(char *flag_byte);
/**
* @brief Get the occupied flag.

0
include/poet/DHT_Types.hpp Normal file → Executable file
View File

116
include/poet/DHT_Wrapper.hpp Normal file → Executable file
View File

@ -1,4 +1,4 @@
// Time-stamp: "Last modified 2023-04-24 16:23:42 mluebke"
// Time-stamp: "Last modified 2023-08-01 13:48:34 mluebke"
/*
** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
@ -23,10 +23,18 @@
#ifndef DHT_WRAPPER_H
#define DHT_WRAPPER_H
#include "DataStructures.hpp"
#include "LookupKey.hpp"
#include "poet/DHT_Types.hpp"
#include "poet/HashFunctions.hpp"
#include "poet/LookupKey.hpp"
#include "poet/Rounding.hpp"
#include <array>
#include <cstdint>
#include <limits>
#include <string>
#include <unordered_map>
#include <utility>
#include <vector>
extern "C" {
@ -37,35 +45,7 @@ extern "C" {
namespace poet {
struct DHT_SCNotation {
std::int8_t exp : 8;
std::int64_t significant : 56;
};
union DHT_Keyelement {
double fp_elemet;
DHT_SCNotation sc_notation;
};
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;
};
/**
* @brief Return user-defined md5sum
*
* This function will calculate a hashsum with the help of md5sum. Therefore the
* md5sum for a given key is calculated and divided into two 64-bit parts. These
* will be XORed and returned as the hash.
*
* @param key_size Size of key in bytes
* @param key Pointer to key
* @return uint64_t Hashsum as an unsigned 64-bit integer
*/
static uint64_t get_md5(int key_size, void *key);
using DHT_Location = std::pair<std::uint32_t, std::uint32_t>;
/**
* @brief C++-Wrapper around DHT implementation
@ -76,6 +56,20 @@ static uint64_t get_md5(int key_size, void *key);
*/
class DHT_Wrapper {
public:
using DHT_ResultObject = struct DHTResobj {
uint32_t length;
std::vector<LookupKey> keys;
std::vector<std::vector<double>> results;
std::vector<double> old_values;
std::vector<bool> needPhreeqc;
std::vector<DHT_Location> locations;
};
static constexpr std::int32_t DHT_KEY_INPUT_CUSTOM =
std::numeric_limits<std::int32_t>::min();
static constexpr int DHT_KEY_SIGNIF_DEFAULT = 5;
/**
* @brief Construct a new dht wrapper object
*
@ -91,8 +85,9 @@ public:
* for key creation.
* @param data_count Count of data entries
*/
DHT_Wrapper(MPI_Comm dht_comm, uint32_t dht_size,
const std::vector<std::uint32_t> &key_indices,
DHT_Wrapper(MPI_Comm dht_comm, std::uint64_t dht_size,
const NamedVector<std::uint32_t> &key_species,
const std::vector<std::int32_t> &key_indices,
uint32_t data_count);
/**
* @brief Destroy the dht wrapper object
@ -105,6 +100,9 @@ public:
*/
~DHT_Wrapper();
DHT_Wrapper &operator=(const DHT_Wrapper &) = delete;
DHT_Wrapper(const DHT_Wrapper &) = delete;
/**
* @brief Check if values of workpackage are stored in DHT
*
@ -123,7 +121,7 @@ public:
*/
auto checkDHT(int length, double dt, const std::vector<double> &work_package,
std::vector<std::uint32_t> &curr_mapping)
-> const poet::DHT_ResultObject &;
-> const DHT_ResultObject &;
/**
* @brief Write simulated values into DHT
@ -183,13 +181,6 @@ public:
*/
auto getHits() { return this->dht_hits; };
/**
* @brief Get the Misses object
*
* @return uint64_t Count of read misses
*/
auto getMisses() { return this->dht_miss; };
/**
* @brief Get the Evictions object
*
@ -197,31 +188,56 @@ public:
*/
auto getEvictions() { return this->dht_evictions; };
void resetCounter() {
this->dht_hits = 0;
this->dht_evictions = 0;
}
void SetSignifVector(std::vector<uint32_t> signif_vec);
void setBaseTotals(const std::array<double, 2> &bt) {
this->base_totals = bt;
auto getDataCount() { return this->data_count; }
auto getCommunicator() { return this->communicator; }
DHT *getDHT() { return this->dht_object; };
DHT_ResultObject &getDHTResults() { return this->dht_results; }
const auto &getKeyElements() const { return this->input_key_elements; }
const auto &getKeySpecies() const { return this->key_species; }
void setBaseTotals(double tot_h, double tot_o) {
this->base_totals = {tot_h, tot_o};
}
std::uint32_t getInputCount() const {
return this->input_key_elements.size();
}
std::uint32_t getOutputCount() const { return this->data_count; }
private:
uint32_t key_count;
uint32_t data_count;
DHT *dht_object;
MPI_Comm communicator;
std::vector<DHT_Keyelement> fuzzForDHT(int var_count, void *key, double dt);
LookupKey fuzzForDHT(int var_count, void *key, double dt);
std::vector<double>
outputToInputAndRates(const std::vector<double> &old_results,
const std::vector<double> &new_results);
std::vector<double>
inputAndRatesToOutput(const std::vector<double> &dht_data);
uint32_t dht_hits = 0;
uint32_t dht_miss = 0;
uint32_t dht_evictions = 0;
std::vector<uint32_t> dht_signif_vector;
std::vector<std::uint32_t> dht_prop_type_vector;
std::vector<std::uint32_t> input_key_elements;
NamedVector<std::uint32_t> key_species;
static constexpr int DHT_KEY_SIGNIF_DEFAULT = 5;
static constexpr int DHT_KEY_SIGNIF_TOTALS = 10;
static constexpr int DHT_KEY_SIGNIF_CHARGE = 3;
std::vector<std::uint32_t> dht_signif_vector;
std::vector<std::uint32_t> dht_prop_type_vector;
std::vector<std::int32_t> input_key_elements;
DHT_ResultObject dht_results;

105
include/poet/DaosKeyValue.h Normal file
View File

@ -0,0 +1,105 @@
/**
* @file DaosKeyValue.h
* @author Nico Sauerbrei (nico.sauerbrei@uni-potsdam.de)
* @brief API to interact with DAOS
* @version 0.1
* @date 01 Jun 2023
*
* This file implements the communication between POET and the DAOS
* Key-Value Store
*/
#ifndef DAOS_KEY_VALUE_H
#define DAOS_KEY_VALUE_H
#include <mpi.h>
#include <stdint.h>
#include <daos.h>
#define DAOS_SUCCESS 0
#define DAOS_ERROR -1
#define DAOS_MPI_ERROR -2
#define DAOS_READ_MISS -3
#define DHT_STATISTICS 1
/**
* Internal struct to store statistics about read and write accesses and also
* read misses and evictions.
* <b>All values will be resetted to zero after a call of
* DHT_print_statistics().</b>
* Internal use only!
*
*/
typedef struct
{
/** Count of writes to specific process this process did. */
int *writes_local;
/** Writes after last call of DHT_print_statistics. */
int old_writes;
/** How many read misses occur? */
int read_misses;
/** How many read hits occur? */
int read_hits;
/** How many buckets where evicted? */
int evictions;
/** How many calls of DHT_write() did this process? */
int w_access;
/** How many calls of DHT_read() did this process? */
int r_access;
} DAOSKV_stats;
/**
* Struct which serves as a handler or so called \a DHT-object. Will
* be created by DHT_create and must be passed as a parameter to every following
* function. Stores all relevant data.
* Do not touch outside DHT functions!
*/
typedef struct
{
/** MPI communicator of all participating processes. */
MPI_Comm communicator;
/** Size of the MPI communicator respectively all participating processes. */
int comm_size;
/** Rank of the process in the MPI communicator. */
int rank;
/** Count of read misses over all time. */
int read_misses;
/** Count of evictions over all time. */
int evictions;
/** Label of the DAOS container.*/
char *cont_label;
/** DAOS pool handle.*/
daos_handle_t poh;
/** DAOS container handle.*/
daos_handle_t coh;
/** DAOS object handle.*/
daos_handle_t oh;
#ifdef DHT_STATISTICS
/** Detailed statistics of the usage of the DHT. */
DAOSKV_stats *stats;
#endif
} DAOSKV;
#ifdef __cplusplus
extern "C"
{
#endif
extern DAOSKV *DAOSKV_create(MPI_Comm comm);
extern int DAOSKV_free(DAOSKV *object);
extern int DAOSKV_write(DAOSKV *object, void *key, int key_size, void *send_data, int send_size);
extern int DAOSKV_read(DAOSKV *object, void *key, int key_size, void *recv_data, int recv_size);
extern int DAOSKV_remove(DAOSKV *object, void *key, int key_size);
extern int DAOSKV_print_statistics(DAOSKV *object);
extern int enumerate_key(DAOSKV *object, int *total_nr, int key_size);
extern struct daos_space get_pool_size(DAOSKV *object);
extern int trim_Space(DAOSKV *object, float deletePercentage, int dataSize, int keySize);
#ifdef __cplusplus
}
#endif
#endif /* DAOS_KEY_VALUE_H */

35
include/poet/DataStructures.hpp Executable file
View File

@ -0,0 +1,35 @@
#ifndef DATASTRUCTURES_H_
#define DATASTRUCTURES_H_
#include <cstddef>
#include <string>
#include <utility>
#include <vector>
namespace poet {
template <typename value_type> class NamedVector {
public:
void insert(std::pair<std::string, value_type> to_insert) {
this->names.push_back(to_insert.first);
this->values.push_back(to_insert.second);
container_size++;
}
bool empty() const { return this->container_size == 0; }
std::size_t size() const { return this->container_size; }
value_type &operator[](std::size_t i) { return values[i]; }
const std::vector<std::string> &getNames() const { return this->names; }
const std::vector<value_type> &getValues() const { return this->values; }
private:
std::size_t container_size{0};
std::vector<std::string> names{};
std::vector<value_type> values{};
};
} // namespace poet
#endif // DATASTRUCTURES_H_

0
include/poet/DiffusionModule.hpp Normal file → Executable file
View File

0
include/poet/Field.hpp Normal file → Executable file
View File

0
include/poet/Grid.hpp Normal file → Executable file
View File

0
include/poet/HashFunctions.hpp Normal file → Executable file
View File

266
include/poet/Interpolation.hpp Executable file
View File

@ -0,0 +1,266 @@
// Time-stamp: "Last modified 2023-08-01 18:10:48 mluebke"
#ifndef INTERPOLATION_H_
#define INTERPOLATION_H_
#include "DHT.h"
#include "DHT_Wrapper.hpp"
#include "DataStructures.hpp"
#include "LookupKey.hpp"
#include "poet/DHT_Wrapper.hpp"
#include "poet/Rounding.hpp"
#include <cassert>
#include <iostream>
#include <list>
#include <memory>
#include <mpi.h>
#include <string>
#include <utility>
extern "C" {
#include "poet/DHT.h"
}
#include "poet/LookupKey.hpp"
#include <cstdint>
#include <functional>
#include <unordered_map>
#include <vector>
namespace poet {
class ProximityHashTable {
public:
using bucket_indicator = std::uint64_t;
ProximityHashTable(uint32_t key_size, uint32_t data_size,
uint32_t entry_count, uint32_t size_per_process,
MPI_Comm communicator);
~ProximityHashTable();
// delete assign and copy operator
ProximityHashTable &operator=(const ProximityHashTable *) = delete;
ProximityHashTable(const ProximityHashTable &) = delete;
struct PHT_Result {
std::uint32_t size;
std::vector<std::vector<double>> in_values;
std::vector<std::vector<double>> out_values;
std::uint32_t getSize() const {
std::uint32_t sum{0};
for (const auto &results : out_values) {
sum += results.size() * sizeof(double);
}
return sum;
}
};
void setSourceDHT(DHT *src) {
this->source_dht = src;
this->dht_key_count = src->key_size / sizeof(Lookup_Keyelement);
this->dht_data_count = src->data_size / sizeof(double);
this->dht_buffer.resize(src->data_size + src->key_size);
}
void writeLocationToPHT(LookupKey key, DHT_Location location);
const PHT_Result &query(const LookupKey &key,
const std::vector<std::uint32_t> &signif,
std::uint32_t min_entries_needed,
std::uint32_t input_count,
std::uint32_t output_count);
std::uint64_t getLocations(const LookupKey &key);
void getEntriesFromLocation(const PHT_Result &locations,
const std::vector<uint32_t> &signif);
void writeStats() { DHT_print_statistics(this->prox_ht); }
DHT *getDHTObject() { return this->prox_ht; }
auto getPHTWriteTime() const -> double { return this->pht_write_t; };
auto getPHTReadTime() const -> double { return this->pht_read_t; };
auto getDHTGatherTime() const -> double { return this->pht_gather_dht_t; };
auto getLocalCacheHits() const -> std::vector<std::uint32_t> {
return this->all_cache_hits;
}
void storeAndResetCounter() {
all_cache_hits.push_back(cache_hits);
cache_hits = 0;
}
#if POET_PHT_ADD
void incrementReadCounter(const LookupKey &key);
#endif
private:
enum { INTERP_CB_OK, INTERP_CB_FULL, INTERP_CB_ALREADY_IN };
static int PHT_callback_function(int in_data_size, void *in_data,
int out_data_size, void *out_data);
static std::vector<double> convertKeysFromDHT(Lookup_Keyelement *keys_in,
std::uint32_t key_size);
static bool similarityCheck(const LookupKey &fine, const LookupKey &coarse,
const std::vector<uint32_t> &signif);
char *bucket_store;
class Cache
: private std::unordered_map<LookupKey, PHT_Result, LookupKeyHasher> {
public:
void operator()(const LookupKey &key, const PHT_Result val);
std::pair<bool, PHT_Result> operator[](const LookupKey &key);
void flush() { this->clear(); }
protected:
private:
static constexpr std::int64_t MAX_CACHE_SIZE = 100E6;
std::int64_t free_mem{MAX_CACHE_SIZE};
std::list<LookupKey> lru_queue;
using lru_iterator = typename std::list<LookupKey>::iterator;
std::unordered_map<LookupKey, lru_iterator, LookupKeyHasher> keyfinder;
};
Cache localCache;
DHT *prox_ht;
std::uint32_t dht_evictions = 0;
DHT *source_dht = nullptr;
PHT_Result lookup_results;
std::vector<char> dht_buffer;
std::uint32_t dht_key_count;
std::uint32_t dht_data_count;
MPI_Comm communicator;
double pht_write_t = 0.;
double pht_read_t = 0.;
double pht_gather_dht_t = 0.;
std::uint32_t cache_hits{0};
std::vector<std::uint32_t> all_cache_hits{};
};
class InterpolationModule {
public:
using InterpFunction = std::vector<double> (*)(
const std::vector<std::int32_t> &, const std::vector<double> &,
const std::vector<std::vector<double>> &,
const std::vector<std::vector<double>> &);
InterpolationModule(std::uint32_t entries_per_bucket,
std::uint64_t size_per_process,
std::uint32_t min_entries_needed, DHT_Wrapper &dht,
const NamedVector<std::uint32_t> &interp_key_signifs,
const std::vector<std::int32_t> &dht_key_indices);
enum result_status { RES_OK, INSUFFICIENT_DATA, NOT_NEEDED };
struct InterpolationResult {
std::vector<std::vector<double>> results;
std::vector<result_status> status;
void ResultsToWP(std::vector<double> &currWP);
};
void setInterpolationFunction(InterpFunction func) {
this->f_interpolate = func;
}
void setMinEntriesNeeded(std::uint32_t entries) {
this->min_entries_needed = entries;
}
auto getMinEntriesNeeded() { return this->min_entries_needed; }
void writePairs(const DHT_Wrapper::DHT_ResultObject &in);
void tryInterpolation(DHT_Wrapper::DHT_ResultObject &dht_results,
std::vector<std::uint32_t> &curr_mapping);
void resultsToWP(std::vector<double> &work_package) const;
auto getPHTWriteTime() const { return pht->getPHTWriteTime(); };
auto getPHTReadTime() const { return pht->getPHTReadTime(); };
auto getDHTGatherTime() const { return pht->getDHTGatherTime(); };
auto getInterpolationTime() const { return this->interp_t; };
auto getInterpolationCount() const -> std::uint32_t {
return this->interpolations;
}
auto getPHTLocalCacheHits() const -> std::vector<std::uint32_t> {
return this->pht->getLocalCacheHits();
}
void resetCounter() {
this->interpolations = 0;
this->pht->storeAndResetCounter();
}
void writePHTStats() { this->pht->writeStats(); }
void dumpPHTState(const std::string &filename) {
DHT_to_file(this->pht->getDHTObject(), filename.c_str());
}
static constexpr std::uint32_t COARSE_DIFF = 2;
static constexpr std::uint32_t COARSE_SIGNIF_DEFAULT =
DHT_Wrapper::DHT_KEY_SIGNIF_DEFAULT - COARSE_DIFF;
private:
void initPHT(std::uint32_t key_count, std::uint32_t entries_per_bucket,
std::uint32_t size_per_process, MPI_Comm communicator);
static std::vector<double> dummy(const std::vector<std::int32_t> &,
const std::vector<double> &,
const std::vector<std::vector<double>> &,
const std::vector<std::vector<double>> &) {
return {};
}
double interp_t = 0.;
std::uint32_t interpolations{0};
InterpFunction f_interpolate = dummy;
std::uint32_t min_entries_needed = 5;
std::unique_ptr<ProximityHashTable> pht;
DHT_Wrapper &dht_instance;
NamedVector<std::uint32_t> key_signifs;
std::vector<std::int32_t> key_indices;
InterpolationResult interp_result;
PHT_Rounder rounder;
LookupKey roundKey(const LookupKey &in_key) {
LookupKey out_key;
for (std::uint32_t i = 0; i < key_indices.size(); i++) {
out_key.push_back(rounder.round(in_key[key_indices[i]], key_signifs[i]));
}
// timestep
out_key.push_back(in_key.back());
return out_key;
}
};
} // namespace poet
#endif // INTERPOLATION_H_

50
include/poet/LookupKey.hpp Executable file
View File

@ -0,0 +1,50 @@
// Time-stamp: "Last modified 2023-07-26 10:20:10 mluebke"
#ifndef LOOKUPKEY_H_
#define LOOKUPKEY_H_
#include "poet/HashFunctions.hpp"
#include <cstdint>
#include <cstring>
#include <vector>
namespace poet {
struct Lookup_SC_notation {
std::int8_t exp : 8;
std::int64_t significant : 56;
};
union Lookup_Keyelement {
double fp_element;
Lookup_SC_notation sc_notation;
bool operator==(const Lookup_Keyelement &other) const {
return std::memcmp(this, &other, sizeof(Lookup_Keyelement)) == 0 ? true
: false;
}
};
class LookupKey : public std::vector<Lookup_Keyelement> {
public:
explicit LookupKey(size_type count);
using std::vector<Lookup_Keyelement>::vector;
std::vector<double> to_double() const;
static Lookup_SC_notation round_from_double(const double in,
std::uint32_t signif);
static double to_double(const Lookup_SC_notation in);
};
struct LookupKeyHasher {
std::uint64_t operator()(const LookupKey &k) const {
const uint32_t key_size = k.size() * sizeof(Lookup_Keyelement);
return poet::Murmur2_64A(key_size, k.data());
}
};
} // namespace poet
#endif // LOOKUPKEY_H_

40
include/poet/RInsidePOET.hpp Normal file → Executable file
View File

@ -1,7 +1,15 @@
#ifndef RPOET_H_
#define RPOET_H_
#include "DataStructures.hpp"
#include <RInside.h>
#include <Rcpp.h>
#include <cstddef>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
class RInsidePOET : public RInside {
public:
@ -14,6 +22,38 @@ public:
RInsidePOET(RInsidePOET const &) = delete;
void operator=(RInsidePOET const &) = delete;
inline bool checkIfExists(const std::string &R_name,
const std::string &where) {
return Rcpp::as<bool>(
this->parseEval("'" + R_name + "' %in% names(" + where + ")"));
}
template <class T> inline T getSTL(const std::string &R_name) {
return Rcpp::as<T>(RInside::operator[](R_name));
}
template <typename value_type>
poet::NamedVector<value_type> wrapNamedVector(const std::string &R_name) {
std::size_t name_size = this->parseEval("length(names(" + R_name + "))");
if (name_size == 0) {
throw std::runtime_error("wrapNamedVector: " + R_name +
" is not a named vector!");
}
auto names = Rcpp::as<std::vector<std::string>>(
this->parseEval("names(" + R_name + ")"));
auto values = Rcpp::as<Rcpp::NumericVector>(this->parseEval(R_name));
poet::NamedVector<value_type> ret_map;
for (std::size_t i = 0; i < names.size(); i++) {
ret_map.insert(
std::make_pair(names[i], static_cast<value_type>(values[i])));
}
return ret_map;
}
private:
RInsidePOET() : RInside(){};
};

92
include/poet/Rounding.hpp Executable file
View File

@ -0,0 +1,92 @@
#ifndef ROUNDING_H_
#define ROUNDING_H_
#include "LookupKey.hpp"
#include <cmath>
#include <cstdint>
namespace poet {
constexpr std::int8_t AQUEOUS_EXP = -13;
template <typename Input, typename Output, typename ConvertTo = double>
class IRounding {
public:
virtual Output round(const Input &, std::uint32_t signif) = 0;
virtual ConvertTo convert(const Output &) = 0;
};
class DHT_Rounder {
public:
Lookup_Keyelement round(const double &value, std::uint32_t signif,
bool is_ho) {
std::int8_t exp =
static_cast<std::int8_t>(std::floor(std::log10(std::fabs(value))));
if (!is_ho) {
if (exp < AQUEOUS_EXP) {
return {.sc_notation = {0, 0}};
}
std::int8_t diff = exp - signif + 1;
if (diff < AQUEOUS_EXP) {
signif -= AQUEOUS_EXP - diff;
}
}
Lookup_Keyelement elem;
elem.sc_notation.exp = exp;
elem.sc_notation.significant =
static_cast<std::int64_t>(value * std::pow(10, signif - exp - 1));
return elem;
}
double convert(const Lookup_Keyelement &key_elem) {
std::int32_t normalized_exp = static_cast<std::int32_t>(
-std::log10(std::fabs(key_elem.sc_notation.significant)));
// add stored exponent to normalized exponent
normalized_exp += key_elem.sc_notation.exp;
// return significant times 10 to the power of exponent
return key_elem.sc_notation.significant * std::pow(10., normalized_exp);
}
};
class PHT_Rounder : public IRounding<Lookup_Keyelement, Lookup_Keyelement> {
public:
Lookup_Keyelement round(const Lookup_Keyelement &value,
std::uint32_t signif) {
Lookup_Keyelement new_val = value;
std::uint32_t diff_signif =
static_cast<std::uint32_t>(
std::ceil(std::log10(std::abs(value.sc_notation.significant)))) -
signif;
new_val.sc_notation.significant = static_cast<int64_t>(
value.sc_notation.significant / std::pow(10., diff_signif));
if (new_val.sc_notation.significant == 0) {
new_val.sc_notation.exp = 0;
}
return new_val;
}
double convert(const Lookup_Keyelement &key_elem) {
std::int32_t normalized_exp = static_cast<std::int32_t>(
-std::log10(std::fabs(key_elem.sc_notation.significant)));
// add stored exponent to normalized exponent
normalized_exp += key_elem.sc_notation.exp;
// return significant times 10 to the power of exponent
return key_elem.sc_notation.significant * std::pow(10., normalized_exp);
}
};
} // namespace poet
#endif // ROUNDING_H_

64
include/poet/SimParams.hpp Normal file → Executable file
View File

@ -24,17 +24,20 @@
#include <cstdint>
#include <string>
#include <string_view>
#include <unordered_map>
#include <vector>
#include "DataStructures.hpp"
#include "RInsidePOET.hpp"
#include "argh.hpp" // Argument handler https://github.com/adishavit/argh
#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 = 1E3;
constexpr uint32_t DHT_SIZE_PER_PROCESS_MB = 1.5E3;
/** Standard work package size */
#define WORK_PACKAGE_SIZE_DEFAULT 5
#define WORK_PACKAGE_SIZE_DEFAULT 32
namespace poet {
@ -70,6 +73,8 @@ typedef struct {
/** 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 {
@ -101,13 +106,24 @@ using DiffusionParams = struct s_DiffusionParams {
s_DiffusionParams(RInside &R);
};
using ChemistryParams = struct s_ChemistryParams {
struct ChemistryParams {
std::string database_path;
std::string input_script;
std::vector<std::string> dht_species;
std::vector<std::uint32_t> dht_signif;
s_ChemistryParams(RInside &R);
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;
void initFromR(RInsidePOET &R);
};
/**
@ -159,7 +175,7 @@ public:
* @return int Returns with 0 if no error occured, otherwise value less than 0
* is returned.
*/
int parseFromCmdl(char *argv[], RInside &R);
int parseFromCmdl(char *argv[], RInsidePOET &R);
/**
* @brief Init std::vector values
@ -193,15 +209,10 @@ public:
*/
auto getDHTSignifVector() const { return this->dht_signif_vector; };
/**
* @brief Return name of DHT snapshot.
*
* Name of the DHT file which is used to initialize the DHT with a previously
* written snapshot.
*
* @return std::string Absolute paht to the DHT snapshot
*/
auto getDHTFile() const { return this->dht_file; };
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
@ -223,22 +234,31 @@ public:
*/
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", "dht-nolog", "P",
"progress"};
const std::set<std::string> paramlist{"work-package-size", "dht-signif",
"dht-strategy", "dht-size",
"dht-snaps", "dht-file"};
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 dht_file;
std::string filesim;
std::string out_dir;
ChemistryParams chem_params;
};
} // namespace poet
#endif // PARSER_H

0
include/poet/argh.hpp Normal file → Executable file
View File

View File

@ -0,0 +1,147 @@
## Time-stamp: "Last modified 2023-08-02 13:59:22 mluebke"
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.217E-09,
"Ba" = 1.E-10,
"Cl" = 2.E-10,
"S" = 6.205E-4,
"Sr" = 6.205E-4,
"Barite" = 0.001,
"Celestite" = 1
)
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.0124,
"H" = 0.00000028904,
#"O" = 55.5087,
"O" = 0.000000165205,
#"Charge" = -1.217E-09,
"Charge" = -3.337E-08,
"Ba" = 1.E-10,
"Cl" = 1.E-10,
"S(6)" = 6.205E-4,
"Sr" = 6.205E-4
)
injection_diff <- list(
list(
#"H" = 111.0124,
"H" = 0.0000002890408,
#"O" = 55.50622,
"O" = 0.00002014464,
#"Charge" = -3.337E-08,
"Charge" = -3.337000004885E-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" = rep(1, n),
## "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, 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" = 10,
"O" = 10,
"Charge" = 3,
"Ba" = 5,
"Cl" = 5,
"S(6)" = 5,
"Sr" = 5
)
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)
)

View File

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

View File

@ -0,0 +1,151 @@
## Time-stamp: "Last modified 2023-07-21 15:04:49 mluebke"
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.217E-09,
"Ba" = 1.E-10,
"Cl" = 2.E-10,
"S" = 6.205E-4,
"Sr" = 6.205E-4,
"Barite" = 0.001,
"Celestite" = 1
)
grid <- list(
n_cells = c(n, m),
s_cells = c(n / 10, m / 10),
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,
"H" = 0.00000028904,
# "O" = 55.5087,
"O" = 0.000000165205,
# "Charge" = -1.217E-09,
"Charge" = -3.337E-08,
"Ba" = 1.E-10,
"Cl" = 1.E-10,
"S(6)" = 6.205E-4,
"Sr" = 6.205E-4
)
injection_diff <- list(
list(
# "H" = 111.0124,
"H" = 0.0000002890408,
# "O" = 55.50622,
"O" = 0.00002014464,
# "Charge" = -3.337E-08,
"Charge" = -3.337000004885E-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, floor(n / 2), floor(m / 2))
## l2 = c(2,80,80),
## l3 = c(2,60,80)
)
boundary <- list(
# "N" = rep(1, n),
"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, 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" = 10,
"O" = 10,
"Charge" = 3,
"Ba" = 5,
"Cl" = 5,
"S(6)" = 5,
"Sr" = 5
)
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)
)

View File

@ -0,0 +1,195 @@
DATABASE
###########################
SOLUTION_MASTER_SPECIES
H H+ -1 H 1.008 # phreeqc/
H(0) H2 0 H # phreeqc/
H(1) H+ -1 0.0 # phreeqc/
E e- 0 0.0 0.0 # phreeqc/
O H2O 0 O 16.0 # phreeqc/
O(0) O2 0 O # phreeqc/
O(-2) H2O 0 0.0 # phreeqc/
Na Na+ 0 Na 22.9898 # phreeqc/
Ba Ba+2 0 Ba 137.34 # phreeqc/
Sr Sr+2 0 Sr 87.62 # phreeqc/
Cl Cl- 0 Cl 35.453 # phreeqc/
S SO4-2 0 SO4 32.064 # phreeqc/
S(6) SO4-2 0 SO4 # phreeqc/
S(-2) HS- 1 S # phreeqc/
SOLUTION_SPECIES
H+ = H+
-gamma 9 0
-dw 9.31e-09
# source: phreeqc
e- = e-
# source: phreeqc
H2O = H2O
# source: phreeqc
Na+ = Na+
-gamma 4.08 0.082
-dw 1.33e-09
-Vm 2.28 -4.38 -4.1 -0.586 0.09 4 0.3 52 -0.00333 0.566
# source: phreeqc
Ba+2 = Ba+2
-gamma 4 0.153
-dw 8.48e-10
-Vm 2.063 -10.06 1.9534 -2.36 0.4218 5 1.58 -12.03 -0.00835 1
# source: phreeqc
Sr+2 = Sr+2
-gamma 5.26 0.121
-dw 7.94e-10
-Vm -0.0157 -10.15 10.18 -2.36 0.86 5.26 0.859 -27 -0.0041 1.97
# source: phreeqc
Cl- = Cl-
-gamma 3.63 0.017
-dw 2.03e-09
-Vm 4.465 4.801 4.325 -2.847 1.748 0 -0.331 20.16 0 1
# source: phreeqc
SO4-2 = SO4-2
-gamma 5 -0.04
-dw 1.07e-09
-Vm 8 2.3 -46.04 6.245 3.82 0 0 0 0 1
# source: phreeqc
H2O = OH- + H+
-analytical_expression 293.29227 0.1360833 -10576.913 -123.73158 0 -6.996455e-05
-gamma 3.5 0
-dw 5.27e-09
-Vm -9.66 28.5 80 -22.9 1.89 0 1.09 0 0 1
# source: phreeqc
2 H2O = O2 + 4 H+ + 4 e-
-log_k -86.08
-delta_h 134.79 kcal
-dw 2.35e-09
-Vm 5.7889 6.3536 3.2528 -3.0417 -0.3943
# source: phreeqc
2 H+ + 2 e- = H2
-log_k -3.15
-delta_h -1.759 kcal
-dw 5.13e-09
-Vm 6.52 0.78 0.12
# source: phreeqc
SO4-2 + H+ = HSO4-
-log_k 1.988
-delta_h 3.85 kcal
-analytical_expression -56.889 0.006473 2307.9 19.8858
-dw 1.33e-09
-Vm 8.2 9.259 2.1108 -3.1618 1.1748 0 -0.3 15 0 1
# source: phreeqc
HS- = S-2 + H+
-log_k -12.918
-delta_h 12.1 kcal
-gamma 5 0
-dw 7.31e-10
# source: phreeqc
SO4-2 + 9 H+ + 8 e- = HS- + 4 H2O
-log_k 33.65
-delta_h -60.140 kcal
-gamma 3.5 0
-dw 1.73e-09
-Vm 5.0119 4.9799 3.4765 -2.9849 1.441
# source: phreeqc
HS- + H+ = H2S
-log_k 6.994
-delta_h -5.30 kcal
-analytical_expression -11.17 0.02386 3279
-dw 2.1e-09
-Vm 7.81 2.96 -0.46
# source: phreeqc
Na+ + OH- = NaOH
-log_k -10
# source: phreeqc
Na+ + SO4-2 = NaSO4-
-log_k 0.7
-delta_h 1.120 kcal
-gamma 5.4 0
-dw 1.33e-09
-Vm 1e-05 16.4 -0.0678 -1.05 4.14 0 6.86 0 0.0242 0.53
# source: phreeqc
Ba+2 + H2O = BaOH+ + H+
-log_k -13.47
-gamma 5 0
# source: phreeqc
Ba+2 + SO4-2 = BaSO4
-log_k 2.7
# source: phreeqc
Sr+2 + H2O = SrOH+ + H+
-log_k -13.29
-gamma 5 0
# source: phreeqc
Sr+2 + SO4-2 = SrSO4
-log_k 2.29
-delta_h 2.08 kcal
-Vm 6.791 -0.9666 6.13 -2.739 -0.001
# source: phreeqc
PHASES
Barite
BaSO4 = Ba+2 + SO4-2
-log_k -9.97
-delta_h 6.35 kcal
-analytical_expression -282.43 -0.08972 5822 113.08
-Vm 52.9
# source: phreeqc
# comment:
Celestite
SrSO4 = Sr+2 + SO4-2
-log_k -6.63
-delta_h -4.037 kcal
-analytical_expression -7.14 0.00611 75 0 0 -1.79e-05
-Vm 46.4
# source: phreeqc
# comment:
RATES
Celestite # Palandri & Kharaka 2004<--------------------------------change me
# PARM(1): reactive surface area
# am: acid mechanism, nm: neutral mechanism, bm: base mechanism
-start
10 sr_i = SR("Celestite") # saturation ratio, (-)<----------change me
20 moles = 0 # init target variable, (mol)
30 IF ((M <= 0) AND (sr_i < 1)) OR (sr_i = 1.0) THEN GOTO 310
40 sa = PARM(1) # reactive surface area, (m2)
100 r = 8.314462 # gas constant, (J K-1 mol-1)
110 dTi = (1 / TK) - (1 / 298.15) # (K-1)
120 ea_am = 23800 # activation energy am, (J mol-1)<-----------change me
130 ea_nm = 0 # activation energy nm, (J mol-1)<-----------change me
140 ea_bm = 0 # activation energy bm, (J mol-1)<-----------change me
150 log_k_am = -5.66 # reaction constant am<-------------------change me
rem log_k_nm = -99 # reaction constant nm<-------------------change me
rem log_k_bm = -99 # reaction constant bm<-------------------change me
180 n_am = 0.109 # H+ reaction order am<-----------------------change me
rem n_bm = 0 # H+ reaction order bm<-----------------------change me
200 am = (10 ^ log_k_am) * EXP(-ea_am * dTi / r) * ACT("H+") ^ n_am
rem nm = (10 ^ log_k_nm) * EXP(-ea_nm * dTi / r)
rem bm = (10 ^ log_k_bm) * EXP(-ea_bm * dTi / r) * ACT("H+") ^ n_bm
300 moles = sa * (am) * (1 - sr_i)
310 save moles * time
-end
Barite # Palandri & Kharaka 2004<-----------------------------------change me
# PARM(1): reactive surface area
# am: acid mechanism, nm: neutral mechanism, bm: base mechanism
-start
10 sr_i = SR("Barite") # saturation ratio, (-)<----------change me
20 moles = 0 # init target variable, (mol)
30 IF ((M <= 0) AND (sr_i < 1)) OR (sr_i = 1.0) THEN GOTO 310
40 sa = PARM(1) # reactive surface area, (m2)
100 r = 8.314462 # gas constant, (J K-1 mol-1)
110 dTi = (1 / TK) - (1 / 298.15) # (K-1)
120 ea_am = 30800 # activation energy am, (J mol-1)<---------change me
130 ea_nm = 30800 # activation energy nm, (J mol-1)<---------change me
rem ea_bm = 0 # activation energy bm, (J mol-1)<---------change me
150 log_k_am = -6.90 # reaction constant am<-----------------change me
160 log_k_nm = -7.90 # reaction constant nm<-----------------change me
rem log_k_bm = -99 # reaction constant bm<-------------------change me
180 n_am = 0.22 # H+ reaction order am<----------------------change me
rem n_bm = 0 # H+ reaction order bm<-----------------------change me
200 am = (10 ^ log_k_am) * EXP(-ea_am * dTi / r) * ACT("H+") ^ n_am
210 nm = (10 ^ log_k_nm) * EXP(-ea_nm * dTi / r)
rem bm = (10 ^ log_k_bm) * EXP(-ea_bm * dTi / r) * ACT("H+") ^ n_bm
300 moles = sa * (am + nm) * (1 - sr_i)
310 save moles * time
-end
END

View File

@ -0,0 +1,205 @@
## Time-stamp: "Last modified 2023-09-05 14:42:20 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, 1, 1)
# 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) ##
#################################################################
## # 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(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
# pht_species = pht_species
)
#################################################################
## Section 4 ##
## Putting all those things together ##
#################################################################
iterations <- 1500
dt <- 500
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))
)

View File

@ -0,0 +1,28 @@
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

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,152 @@
## Time-stamp: "Last modified 2023-08-02 13:59:02 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
)
chemistry <- list(
database = database,
input_script = input_script,
dht_species = dht_species
)
#################################################################
## 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
)

View File

@ -0,0 +1,152 @@
## Time-stamp: "Last modified 2023-08-02 13:59:12 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(2, 1),
type = types[1]
)
##################################################################
## Section 2 ##
## Diffusion parameters and boundary conditions ##
##################################################################
## initial conditions
init_diffu <- list(
"H" = 0.000211313883539788,
"O" = 0.00398302904424952,
"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" = 0.0001540445,
"O" = 0.002148006,
"Charge" = 1.90431e-16,
"C(4)" = 0,
"Ca" = 0,
"Cl" = 0.002,
"Mg" = 0.001
),
list(
"H" = 0.0001610193,
"O" = 0.002386934,
"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
)
chemistry <- list(
database = database,
input_script = input_script,
dht_species = dht_species
)
#################################################################
## 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)
)

View File

@ -0,0 +1,28 @@
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

View File

@ -0,0 +1,167 @@
## Time-stamp: "Last modified 2023-08-02 13:47:06 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
#)
chemistry <- list(
database = database,
input_script = input_script,
dht_species = dht_species
# 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))
)

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,39 @@
## Time-stamp: "Last modified 2023-03-21 11:49:43 mluebke"
KNOBS
-logfile false
-iterations 10000
-convergence_tolerance 1E-12
-step_size 2
-pe_step_size 2
SELECTED_OUTPUT
-reset false
-high_precision true
-solution true
-state true
-step true
-pH true
-pe true
-ionic_strength true
-water true
SOLUTION 1
temp 13
units mol/kgw
pH 7.06355
pe -2.626517
C(4) 0.001990694
Ca 0.02172649
Cl 0.3227673 charge
Fe 0.0001434717
K 0.001902357
Mg 0.01739704
Na 0.2762882
S(6) 0.01652701
Sr 0.0004520361
U(4) 8.147792e-12
U(6) 2.237946e-09
-water 0.00133
EXCHANGE 1
-equil 1
Z 0.0012585
Y 0.0009418
END

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,56 @@
## Time-stamp: "Last modified 2023-02-27 14:31:11 delucia"
KNOBS
-logfile false
-iterations 10000
-convergence_tolerance 1E-12
-step_size 2
-pe_step_size 2
SELECTED_OUTPUT
-reset false
-high_precision true
-solution true
-state true
-step true
-pH true
-pe true
-ionic_strength true
-water true
USER_PUNCH
-head total_o total_h cb C(-4) C(4) Ca Cl Fe(2) Fe(3) H(0) K Mg Na S(-2) S(2) S(4) S(6) Sr U(4) U(5) U(6) UO2(am,hyd) KdU
-start
5 w=TOT("water")
10 PUNCH TOTMOLE("O"), TOTMOLE("H"), CHARGE_BALANCE, w*TOT("C(-4)"), w*TOT("C(4)"), w*TOT("Ca"), w*TOT("Cl"), w*TOT("Fe(2)"), w*TOT("Fe(3)"), w*TOT("H(0)"), w*TOT("K"), w*TOT("Mg"), w*TOT("Na"), w*TOT("S(-2)"), w*TOT("S(2)"), w*TOT("S(4)"), w*TOT("S(6)"), w*TOT("Sr"), w*TOT("U(4)"), w*TOT("U(5)"), w*TOT("U(6)"), EQUI("UO2(am,hyd)")
20 PUNCH ((SURF("U, Ill")+SURF("U, Mll")+SURF("U, Kln")+EDL("U, Ill")+EDL("U, Mll")+EDL("U, Kln"))/((TOT("U")*1.01583)))/(0.002251896406*1000)
-end
SOLUTION 1
temp 13
units mol/kgw
pH 7.06355
pe -2.626517
C(4) 0.001990694
Ca 0.02172649
Cl 0.3227673 charge
Fe 0.0001434717
K 0.001902357
Mg 0.01739704
Na 0.2762882
S(6) 0.01652701
Sr 0.0004520361
U(4) 8.147792e-12
U(6) 2.237946e-09
-water 0.00133
SURFACE 1
-equil 1
-sites_units density
-donnan 4.9e-10
Kln_aOH 1.155 11. 5.0518
Kln_siOH 1.155
Ill_sOH 0.05 100. 5.5931
Ill_wOH 2.26
Mll_sOH 0.05 100. 1.0825
Mll_wOH 2.26
EXCHANGE 1
-equil 1
Z 0.0012585
Y 0.0009418
END

View File

@ -0,0 +1,140 @@
## Time-stamp: "Last modified 2023-08-02 13:59:35 mluebke"
database <- normalizePath("./SMILE_2021_11_01_TH.dat")
input_script <- normalizePath("./ExBase.pqi")
cat(paste(":: R This is a test 1\n"))
#################################################################
## Section 1 ##
## Grid initialization ##
#################################################################
n <- 100
m <- 100
types <- c("scratch", "phreeqc", "rds")
init_cell <- list(H = 1.476571028625e-01,
O = 7.392297218936e-02,
Charge = -1.765225732724e-18,
`C(-4)` = 2.477908970828e-21,
`C(4)` = 2.647623016916e-06,
Ca = 2.889623169138e-05,
Cl = 4.292806181039e-04,
`Fe(2)` =1.908142472666e-07,
`Fe(3)` =3.173306589931e-12,
`H(0)` =2.675642675119e-15,
K = 2.530134809667e-06,
Mg =2.313806319294e-05,
Na =3.674633059628e-04,
`S(-2)` = 8.589766637180e-15,
`S(2)` = 1.205284362720e-19,
`S(4)` = 9.108958772790e-18,
`S(6)` = 2.198092329098e-05,
Sr = 6.012080128154e-07,
`U(4)` = 1.039668623852e-14,
`U(5)` = 1.208394829796e-15,
`U(6)` = 2.976409147150e-12)
grid <- list(
n_cells = c(n, m),
s_cells = c(1, 1),
type = "scratch"
)
##################################################################
## 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
setup <- list(
grid = grid,
diffusion = diffusion,
chemistry = chemistry,
iterations = iterations,
timesteps = rep(dt, iterations),
store_result = TRUE
)

View File

@ -0,0 +1,141 @@
## Time-stamp: "Last modified 2023-08-02 13:59:44 mluebke"
database <- normalizePath("../share/poet/bench/surfex/SMILE_2021_11_01_TH.dat")
input_script <- normalizePath("../share/poet/bench/surfex/SurfExBase.pqi")
cat(paste(":: R This is a test 1\n"))
#################################################################
## Section 1 ##
## Grid initialization ##
#################################################################
n <- 10
m <- 10
types <- c("scratch", "phreeqc", "rds")
init_cell <- list(H = 1.476571028625e-01,
O = 7.392297218936e-02,
Charge = -1.765225732724e-18,
`C(-4)` = 2.477908970828e-21,
`C(4)` = 2.647623016916e-06,
Ca = 2.889623169138e-05,
Cl = 4.292806181039e-04,
`Fe(2)` =1.908142472666e-07,
`Fe(3)` =3.173306589931e-12,
`H(0)` =2.675642675119e-15,
K = 2.530134809667e-06,
Mg =2.313806319294e-05,
Na =3.674633059628e-04,
`S(-2)` = 8.589766637180e-15,
`S(2)` = 1.205284362720e-19,
`S(4)` = 9.108958772790e-18,
`S(6)` = 2.198092329098e-05,
Sr = 6.012080128154e-07,
`U(4)` = 1.039668623852e-14,
`U(5)` = 1.208394829796e-15,
`U(6)` = 2.976409147150e-12)
grid <- list(
n_cells = c(n, m),
s_cells = c(1, 1),
type = "scratch"
)
##################################################################
## 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
setup <- list(
grid = grid,
diffusion = diffusion,
chemistry = chemistry,
iterations = iterations,
timesteps = rep(dt, iterations),
store_result = TRUE,
out_save = c(5, iterations)
)

50
src/CMakeLists.txt Normal file → Executable file
View File

@ -1,15 +1,55 @@
option(POET_DAOSKV "Build with DAOS, using the KV-API" OFF)
option(POET_DAOSOB "Build with DAOS, using the Object-API" OFF)
file(GLOB_RECURSE poet_lib_SRC
CONFIGURE_DEPENDS
"*.cpp" "*.c")
find_library(MATH_LIBRARY m)
find_library(CRYPTO_LIBRARY crypto)
if(POET_DAOSKV OR POET_DAOSOB)
#list(REMOVE_ITEM poet_lib_SRC "${CMAKE_CURRENT_SOURCE_DIR}/ChemistryModule/SurrogateModels/DHT.c")
list(REMOVE_ITEM poet_lib_SRC "${CMAKE_CURRENT_SOURCE_DIR}/ChemistryModule/SurrogateModels/DHT_Wrapper.cpp")
if (POET_DAOSKV)
message(NOTICE "Using DAOS KV")
list(REMOVE_ITEM poet_lib_SRC "${CMAKE_CURRENT_SOURCE_DIR}/ChemistryModule/DAOS/OB/DaosKeyValue.c")
else ()
message(NOTICE "Using DAOS OB")
list(REMOVE_ITEM poet_lib_SRC "${CMAKE_CURRENT_SOURCE_DIR}/ChemistryModule/DAOS/KV/DaosKeyValue.c")
endif()
# ML: Fuer Nico: Includes + linking fuer Daos bereitstellen
find_library(DAOS_LIB libdaos.so
PATH_SUFFIXES lib lib64
)
find_path(DAOS_INCLUDE daos.h)
add_library(DAOS INTERFACE IMPORTED)
set_target_properties(DAOS PROPERTIES
INTERFACE_LINK_LIBRARIES "${DAOS_LIB}"
INTERFACE_INCLUDE_DIRECTORIES "${DAOS_INCLUDE}"
)
set(PHREEQCRM_BUILD_MPI OFF CACHE BOOL "" FORCE)
else()
list(REMOVE_ITEM poet_lib_SRC "${CMAKE_CURRENT_SOURCE_DIR}/ChemistryModule/DAOS/DHT_Wrapper.cpp")
list(REMOVE_ITEM poet_lib_SRC "${CMAKE_CURRENT_SOURCE_DIR}/ChemistryModule/DAOS/KV/DaosKeyValue.c")
list(REMOVE_ITEM poet_lib_SRC "${CMAKE_CURRENT_SOURCE_DIR}/ChemistryModule/DAOS/OB/DaosKeyValue.c")
endif()
foreach(src IN LISTS poet_lib_SRC)
message(STATUS " ${src}")
endforeach()
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_CXX ${MATH_LIBRARY} RRuntime PhreeqcRM tug)
if(POET_DAOSKV OR POET_DAOSOB)
target_link_libraries(poet_lib PUBLIC DAOS)
endif()
target_compile_definitions(poet_lib PUBLIC STRICT_R_HEADERS OMPI_SKIP_MPICXX)
mark_as_advanced(PHREEQCRM_BUILD_MPI PHREEQCRM_DISABLE_OPENMP)
@ -20,3 +60,9 @@ option(POET_DHT_DEBUG "Build with DHT debug info" OFF)
if(POET_DHT_DEBUG)
target_compile_definitions(poet_lib PRIVATE DHT_STATISTICS)
endif()
option(POET_PHT_ADDITIONAL_INFO "Enables additional information in the PHT" OFF)
if (POET_PHT_ADDITIONAL_INFO)
target_compile_definitions(poet_lib PRIVATE POET_PHT_ADD)
endif()

385
src/ChemistryModule/ChemistryModule.cpp Normal file → Executable file
View File

@ -1,20 +1,167 @@
#include "poet/ChemistryModule.hpp"
#include "PhreeqcRM.h"
#include "poet/DHT_Wrapper.hpp"
#include "poet/Interpolation.hpp"
#include <algorithm>
#include <cassert>
#include <cstddef>
#include <cstdint>
#include <map>
#include <memory>
#include <mpi.h>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
constexpr uint32_t MB_FACTOR = 1E6;
std::vector<double>
inverseDistanceWeighting(const std::vector<std::int32_t> &to_calc,
const std::vector<double> &from,
const std::vector<std::vector<double>> &input,
const std::vector<std::vector<double>> &output) {
std::vector<double> results = from;
const std::uint32_t buffer_size = input.size() + 1;
double buffer[buffer_size];
double from_rescaled;
const std::uint32_t data_set_n = input.size();
double rescaled[to_calc.size()][data_set_n + 1];
double weights[data_set_n];
// rescaling over all key elements
for (int key_comp_i = 0; key_comp_i < to_calc.size(); key_comp_i++) {
const auto output_comp_i = to_calc[key_comp_i];
// rescale input between 0 and 1
for (int point_i = 0; point_i < data_set_n; point_i++) {
rescaled[key_comp_i][point_i] = input[point_i][key_comp_i];
}
rescaled[key_comp_i][data_set_n] = from[output_comp_i];
const double min = *std::min_element(rescaled[key_comp_i],
rescaled[key_comp_i] + data_set_n + 1);
const double max = *std::max_element(rescaled[key_comp_i],
rescaled[key_comp_i] + data_set_n + 1);
for (int point_i = 0; point_i < data_set_n; point_i++) {
rescaled[key_comp_i][point_i] =
((max - min) != 0
? (rescaled[key_comp_i][point_i] - min) / (max - min)
: 0);
}
rescaled[key_comp_i][data_set_n] =
((max - min) != 0 ? (from[output_comp_i] - min) / (max - min) : 0);
}
// calculate distances for each data set
double inv_sum = 0;
for (int point_i = 0; point_i < data_set_n; point_i++) {
double distance = 0;
for (int key_comp_i = 0; key_comp_i < to_calc.size(); key_comp_i++) {
distance += std::pow(
rescaled[key_comp_i][point_i] - rescaled[key_comp_i][data_set_n], 2);
}
weights[point_i] = 1 / std::sqrt(distance);
assert(!std::isnan(weights[point_i]));
inv_sum += weights[point_i];
}
assert(!std::isnan(inv_sum));
// actual interpolation
// bool has_h = false;
// bool has_o = false;
for (int key_comp_i = 0; key_comp_i < to_calc.size(); key_comp_i++) {
const auto output_comp_i = to_calc[key_comp_i];
double key_delta = 0;
// if (interp_i == 0) {
// has_h = true;
// }
// if (interp_i == 1) {
// has_o = true;
// }
for (int j = 0; j < data_set_n; j++) {
key_delta += weights[j] * input[j][key_comp_i];
}
key_delta /= inv_sum;
results[output_comp_i] = from[output_comp_i] + key_delta;
}
// if (!has_h) {
// double new_val = 0;
// for (int j = 0; j < data_set_n; j++) {
// new_val += weights[j] * output[j][0];
// }
// results[0] = new_val / inv_sum;
// }
// if (!has_h) {
// double new_val = 0;
// for (int j = 0; j < data_set_n; j++) {
// new_val += weights[j] * output[j][1];
// }
// results[1] = new_val / inv_sum;
// }
// for (std::uint32_t i = 0; i < to_calc.size(); i++) {
// const std::uint32_t interp_i = to_calc[i];
// // rescale input between 0 and 1
// for (int j = 0; j < input.size(); j++) {
// buffer[j] = input[j].at(i);
// }
// buffer[buffer_size - 1] = from[interp_i];
// const double min = *std::min_element(buffer, buffer + buffer_size);
// const double max = *std::max_element(buffer, buffer + buffer_size);
// for (int j = 0; j < input.size(); j++) {
// buffer[j] = ((max - min) != 0 ? (buffer[j] - min) / (max - min) : 1);
// }
// from_rescaled =
// ((max - min) != 0 ? (from[interp_i] - min) / (max - min) : 0);
// double inv_sum = 0;
// // calculate distances for each point
// for (int i = 0; i < input.size(); i++) {
// const double distance = std::pow(buffer[i] - from_rescaled, 2);
// buffer[i] = distance > 0 ? (1 / std::sqrt(distance)) : 0;
// inv_sum += buffer[i];
// }
// // calculate new values
// double new_val = 0;
// for (int i = 0; i < output.size(); i++) {
// new_val += buffer[i] * output[i][interp_i];
// }
// results[interp_i] = new_val / inv_sum;
// if (std::isnan(results[interp_i])) {
// std::cout << "nan with new_val = " << output[0][i] << std::endl;
// }
// }
return results;
}
poet::ChemistryModule::ChemistryModule(uint32_t nxyz, uint32_t wp_size,
std::uint32_t maxiter,
const ChemistryParams &chem_param,
MPI_Comm communicator)
: PhreeqcRM(nxyz, 1), group_comm(communicator), wp_size(wp_size) {
: PhreeqcRM(nxyz, 1), group_comm(communicator), wp_size(wp_size),
params(chem_param) {
MPI_Comm_size(communicator, &this->comm_size);
MPI_Comm_rank(communicator, &this->comm_rank);
@ -24,7 +171,10 @@ poet::ChemistryModule::ChemistryModule(uint32_t nxyz, uint32_t wp_size,
if (!is_sequential && is_master) {
MPI_Bcast(&wp_size, 1, MPI_UINT32_T, 0, this->group_comm);
MPI_Bcast(&maxiter, 1, MPI_UINT32_T, 0, this->group_comm);
}
this->file_pad = std::ceil(std::log10(maxiter + 1));
}
poet::ChemistryModule::~ChemistryModule() {
@ -34,11 +184,15 @@ poet::ChemistryModule::~ChemistryModule() {
}
poet::ChemistryModule
poet::ChemistryModule::createWorker(MPI_Comm communicator) {
poet::ChemistryModule::createWorker(MPI_Comm communicator,
const ChemistryParams &chem_param) {
uint32_t wp_size;
MPI_Bcast(&wp_size, 1, MPI_UINT32_T, 0, communicator);
return ChemistryModule(wp_size, wp_size, 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) {
@ -122,7 +276,9 @@ void poet::ChemistryModule::initializeField(const Field &trans_field) {
}
// now sort the new values
std::sort(new_solution_names.begin() + 4, new_solution_names.end());
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;
@ -133,8 +289,6 @@ void poet::ChemistryModule::initializeField(const Field &trans_field) {
this->prop_names = std::move(new_prop_names);
this->prop_count = prop_names.size();
this->SetPOETSolutionNames(new_solution_names);
if (is_master) {
this->n_cells = trans_field.GetRequestedVecSize();
chem_field = Field(n_cells);
@ -151,145 +305,186 @@ void poet::ChemistryModule::initializeField(const Field &trans_field) {
std::vector<std::vector<double>> initial_values;
for (int i = 0; i < phreeqc_init.size(); i++) {
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));
}
chem_field.InitFromVec(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::SetDHTEnabled(
bool enable, uint32_t size_mb,
const std::vector<std::string> &key_species) {
void poet::ChemistryModule::initializeDHT(
uint32_t size_mb, const NamedVector<std::uint32_t> &key_species) {
constexpr uint32_t MB_FACTOR = 1E6;
std::vector<std::uint32_t> key_inidices;
this->dht_enabled = true;
MPI_Comm dht_comm;
if (this->is_master) {
int ftype = CHEM_DHT_ENABLE;
PropagateFunctionType(ftype);
ChemBCast(&enable, 1, MPI_CXX_BOOL);
ChemBCast(&size_mb, 1, MPI_UINT32_T);
key_inidices = parseDHTSpeciesVec(key_species);
int vec_size = key_inidices.size();
ChemBCast(&vec_size, 1, MPI_INT);
ChemBCast(key_inidices.data(), key_inidices.size(), MPI_UINT32_T);
} else {
int vec_size;
ChemBCast(&vec_size, 1, MPI_INT);
key_inidices.resize(vec_size);
ChemBCast(key_inidices.data(), vec_size, MPI_UINT32_T);
MPI_Comm_split(this->group_comm, MPI_UNDEFINED, this->comm_rank, &dht_comm);
return;
}
this->dht_enabled = enable;
if (enable) {
MPI_Comm dht_comm;
if (this->is_master) {
MPI_Comm_split(this->group_comm, MPI_UNDEFINED, this->comm_rank,
&dht_comm);
return;
}
if (!this->is_master) {
MPI_Comm_split(this->group_comm, 1, this->comm_rank, &dht_comm);
auto map_copy = key_species;
if (key_species.empty()) {
const auto &all_species = this->prop_names;
for (std::size_t i = 0; i < all_species.size(); i++) {
map_copy.insert(std::make_pair(all_species[i],
DHT_Wrapper::DHT_KEY_SIGNIF_DEFAULT));
}
}
auto key_indices = parseDHTSpeciesVec(key_species, this->prop_names);
if (this->dht) {
delete this->dht;
}
const uint32_t dht_size = size_mb * MB_FACTOR;
const std::uint64_t dht_size = size_mb * MB_FACTOR;
this->dht =
new DHT_Wrapper(dht_comm, dht_size, key_inidices, this->prop_count);
// this->dht->setBaseTotals(this->base_totals);
this->dht = new DHT_Wrapper(dht_comm, dht_size, map_copy, key_indices,
this->prop_count);
this->dht->setBaseTotals(base_totals.at(0), base_totals.at(1));
}
}
inline std::vector<std::uint32_t> poet::ChemistryModule::parseDHTSpeciesVec(
const std::vector<std::string> &species_vec) const {
std::vector<uint32_t> species_indices;
inline std::vector<std::int32_t> poet::ChemistryModule::parseDHTSpeciesVec(
const NamedVector<std::uint32_t> &key_species,
const std::vector<std::string> &to_compare) const {
std::vector<int32_t> species_indices;
species_indices.reserve(key_species.size());
if (species_vec.empty()) {
species_indices.resize(this->prop_count);
int i = 0;
std::generate(species_indices.begin(), species_indices.end(),
[&] { return i++; });
return species_indices;
}
species_indices.reserve(species_vec.size());
for (const auto &name : species_vec) {
auto it = std::find(this->prop_names.begin(), this->prop_names.end(), name);
if (it == prop_names.end()) {
throw std::runtime_error(
"DHT species name was not found in prop name vector!");
for (const auto &species : key_species.getNames()) {
auto it = std::find(to_compare.begin(), to_compare.end(), species);
if (it == to_compare.end()) {
species_indices.push_back(DHT_Wrapper::DHT_KEY_INPUT_CUSTOM);
continue;
}
const std::uint32_t index = it - prop_names.begin();
const std::uint32_t index = it - to_compare.begin();
species_indices.push_back(index);
}
return species_indices;
}
void poet::ChemistryModule::SetDHTSnaps(int type, const std::string &out_dir) {
void poet::ChemistryModule::BCastStringVec(std::vector<std::string> &io) {
if (this->is_master) {
int ftype = CHEM_DHT_SNAPS;
PropagateFunctionType(ftype);
int vec_size = io.size();
ChemBCast(&vec_size, 1, MPI_INT);
int str_size = out_dir.size();
for (const auto &value : io) {
int buf_size = value.size() + 1;
ChemBCast(&buf_size, 1, MPI_INT);
ChemBCast(const_cast<char *>(value.c_str()), buf_size, MPI_CHAR);
}
} else {
int vec_size;
ChemBCast(&vec_size, 1, MPI_INT);
ChemBCast(&type, 1, MPI_INT);
ChemBCast(&str_size, 1, MPI_INT);
ChemBCast(const_cast<char *>(out_dir.data()), str_size, MPI_CHAR);
io.resize(vec_size);
for (int i = 0; i < vec_size; i++) {
int buf_size;
ChemBCast(&buf_size, 1, MPI_INT);
char buf[buf_size];
ChemBCast(buf, buf_size, MPI_CHAR);
io[i] = std::string{buf};
}
}
this->dht_snaps_type = type;
this->dht_file_out_dir = out_dir;
}
void poet::ChemistryModule::SetDHTSignifVector(
std::vector<uint32_t> &signif_vec) {
void poet::ChemistryModule::setDHTSnapshots(int type,
const std::string &out_dir) {
if (this->is_master) {
int ftype = CHEM_DHT_SIGNIF_VEC;
PropagateFunctionType(ftype);
int data_count = signif_vec.size();
ChemBCast(&data_count, 1, MPI_INT);
ChemBCast(signif_vec.data(), signif_vec.size(), MPI_UINT32_T);
return;
}
int data_count;
ChemBCast(&data_count, 1, MPI_INT);
signif_vec.resize(data_count);
ChemBCast(signif_vec.data(), data_count, MPI_UINT32_T);
this->dht->SetSignifVector(signif_vec);
this->dht_file_out_dir = out_dir;
this->dht_snaps_type = type;
}
void poet::ChemistryModule::ReadDHTFile(const std::string &input_file) {
void poet::ChemistryModule::setDHTReadFile(const std::string &input_file) {
if (this->is_master) {
int ftype = CHEM_DHT_READ_FILE;
PropagateFunctionType(ftype);
int str_size = input_file.size();
ChemBCast(&str_size, 1, MPI_INT);
ChemBCast(const_cast<char *>(input_file.data()), str_size, MPI_CHAR);
return;
}
if (!this->dht_enabled) {
throw std::runtime_error("DHT file cannot be read. DHT is not enabled.");
if (!input_file.empty()) {
WorkerReadDHTDump(input_file);
}
}
void poet::ChemistryModule::initializeInterp(
std::uint32_t bucket_size, std::uint32_t size_mb, std::uint32_t min_entries,
const NamedVector<std::uint32_t> &key_species) {
if (!this->is_master) {
WorkerReadDHTDump(input_file);
constexpr uint32_t MB_FACTOR = 1E6;
assert(this->dht);
this->interp_enabled = true;
auto map_copy = key_species;
if (key_species.empty()) {
map_copy = this->dht->getKeySpecies();
for (std::size_t i = 0; i < map_copy.size(); i++) {
const std::uint32_t signif =
map_copy[i] - (map_copy[i] > InterpolationModule::COARSE_DIFF
? InterpolationModule::COARSE_DIFF
: 0);
map_copy[i] = signif;
}
}
auto key_indices =
parseDHTSpeciesVec(key_species, dht->getKeySpecies().getNames());
if (this->interp) {
this->interp.reset();
}
const uint64_t pht_size = size_mb * MB_FACTOR;
interp = std::make_unique<poet::InterpolationModule>(
bucket_size, pht_size, min_entries, *(this->dht), map_copy,
key_indices);
interp->setInterpolationFunction(inverseDistanceWeighting);
}
}

View File

@ -0,0 +1,313 @@
// Time-stamp: "Last modified 2023-08-01 13:41:57 mluebke"
/*
** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
** Potsdam)
**
** Copyright (C) 2018-2021 Marco De Lucia (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 "poet/DHT_Wrapper.hpp"
#include "poet/DHT_Types.hpp"
#include "poet/HashFunctions.hpp"
#include "poet/Interpolation.hpp"
#include "poet/LookupKey.hpp"
#include "poet/Rounding.hpp"
#include "poet/DaosKeyValue.h"
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <cstring>
#include <iostream>
#include <stdexcept>
#include <vector>
using namespace std;
namespace poet {
uint32_t data_size = 0;
uint32_t key_size = 0;
DAOSKV *daosKV_object;
DHT_Wrapper::DHT_Wrapper(MPI_Comm dht_comm, std::uint64_t dht_size,
const NamedVector<std::uint32_t> &key_species,
const std::vector<std::int32_t> &key_indices,
uint32_t data_count)
: key_count(key_indices.size()), data_count(data_count),
input_key_elements(key_indices), communicator(dht_comm),
key_species(key_species) {
// initialize DHT object
// key size = count of key elements + timestep
key_size = (key_count + 1) * sizeof(Lookup_Keyelement);
data_size =
(data_count + input_key_elements.size()) * sizeof(double);
uint32_t buckets_per_process =
static_cast<std::uint32_t>(dht_size / (data_size + key_size));
daosKV_object = DAOSKV_create(dht_comm);
dht_signif_vector = key_species.getValues();
// this->dht_signif_vector.resize(key_size, DHT_KEY_SIGNIF_DEFAULT);
this->dht_prop_type_vector.resize(key_count, DHT_TYPE_DEFAULT);
const auto key_names = key_species.getNames();
auto tot_h = std::find(key_names.begin(), key_names.end(), "H");
if (tot_h != key_names.end()) {
this->dht_prop_type_vector[tot_h - key_names.begin()] = DHT_TYPE_TOTAL;
}
auto tot_o = std::find(key_names.begin(), key_names.end(), "O");
if (tot_o != key_names.end()) {
this->dht_prop_type_vector[tot_o - key_names.begin()] = DHT_TYPE_TOTAL;
}
auto charge = std::find(key_names.begin(), key_names.end(), "Charge");
if (charge != key_names.end()) {
this->dht_prop_type_vector[charge - key_names.begin()] = DHT_TYPE_CHARGE;
}
}
DHT_Wrapper::~DHT_Wrapper() {
// free DHT
DAOSKV_free(daosKV_object);
}
auto DHT_Wrapper::checkDHT(int length, double dt,
const std::vector<double> &work_package,
std::vector<std::uint32_t> &curr_mapping)
-> const DHT_ResultObject & {
dht_results.length = length;
dht_results.keys.resize(length);
dht_results.results.resize(length);
dht_results.needPhreeqc.resize(length);
std::vector<double> bucket_writer(this->data_count +
input_key_elements.size());
std::vector<std::uint32_t> new_mapping;
// loop over every grid cell contained in work package
for (int i = 0; i < length; i++) {
// point to current grid cell
void *key = (void *)&(work_package[i * this->data_count]);
auto &data = dht_results.results[i];
auto &key_vector = dht_results.keys[i];
// 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
int res = DAOSKV_read(daosKV_object,
key_vector.data(),key_size,
bucket_writer.data(), data_size);
switch (res) {
case DAOS_SUCCESS:
dht_results.results[i] = inputAndRatesToOutput(bucket_writer);
dht_results.needPhreeqc[i] = false;
this->dht_hits++;
break;
case DAOS_READ_MISS:
dht_results.needPhreeqc[i] = true;
new_mapping.push_back(curr_mapping[i]);
dht_results.results[i] = std::vector<double>{
&work_package[i * this->data_count],
&work_package[i * this->data_count] + this->data_count};
// HACK: apply normalization to total H and O in results field of DHT
// dht_results.results[i][0] -= base_totals[0];
// dht_results.results[i][1] -= base_totals[1];
break;
}
}
curr_mapping = std::move(new_mapping);
dht_results.old_values = work_package;
return dht_results;
}
void DHT_Wrapper::fillDHT(int length, const std::vector<double> &work_package) {
// loop over every grid cell contained in work package
dht_results.locations.resize(length);
for (int i = 0; i < length; i++) {
// If true grid cell was simulated, needs to be inserted into dht
if (dht_results.needPhreeqc[i]) {
// check if calcite or dolomite is absent and present, resp.n and vice
// versa in input/output. If this is the case -> Do not write to DHT!
// HACK: hardcoded, should be fixed!
if ((dht_results.old_values[i * this->data_count + 7] == 0) !=
(work_package[i * this->data_count + 7] == 0)) {
dht_results.needPhreeqc[i] = false;
continue;
}
if ((dht_results.old_values[i * this->data_count + 9] == 0) !=
(work_package[i * this->data_count + 9] == 0)) {
dht_results.needPhreeqc[i] = false;
continue;
}
uint32_t proc, index;
const auto &key = dht_results.keys[i];
const auto curr_old_data = std::vector<double>(
dht_results.old_values.begin() + (i * this->data_count),
dht_results.old_values.begin() + ((i + 1) * this->data_count));
const auto curr_new_data = std::vector<double>(
work_package.begin() + (i * this->data_count),
work_package.begin() + ((i + 1) * this->data_count));
const auto data = outputToInputAndRates(curr_old_data, curr_new_data);
// void *data = (void *)&(work_package[i * this->data_count]);
// fuzz data (round, logarithm etc.)
// insert simulated data with fuzzed key into DHT
int res = DAOSKV_write(daosKV_object,
(void *) key.data(), key_size,
(void *) data.data(), data_size);
dht_results.locations[i] = {proc, index};
// if data was successfully written ...
if ((res != DHT_SUCCESS) && (res == DHT_WRITE_SUCCESS_WITH_EVICTION)) {
dht_evictions++;
}
}
}
}
std::vector<double>
DHT_Wrapper::outputToInputAndRates(const std::vector<double> &old_results,
const std::vector<double> &new_results) {
const int prefix_size = this->input_key_elements.size();
std::vector<double> output(prefix_size + this->data_count);
std::copy(new_results.begin(), new_results.end(),
output.begin() + prefix_size);
for (int i = 0; i < prefix_size; i++) {
const int data_elem_i = input_key_elements[i];
output[i] = old_results[data_elem_i];
output[prefix_size + data_elem_i] -= old_results[data_elem_i];
}
return output;
}
std::vector<double>
DHT_Wrapper::inputAndRatesToOutput(const std::vector<double> &dht_data) {
const int prefix_size = this->input_key_elements.size();
std::vector<double> output{dht_data.begin() + prefix_size, dht_data.end()};
for (int i = 0; i < prefix_size; i++) {
const int data_elem_i = input_key_elements[i];
output[data_elem_i] += dht_data[i];
}
return output;
}
void DHT_Wrapper::resultsToWP(std::vector<double> &work_package) {
for (int i = 0; i < dht_results.length; i++) {
if (!dht_results.needPhreeqc[i]) {
std::copy(dht_results.results[i].begin(), dht_results.results[i].end(),
work_package.begin() + (data_count * i));
}
}
}
int DHT_Wrapper::tableToFile(const char *filename) {
int res = 0; //DHT_to_file(dht_object, filename);
return res;
}
int DHT_Wrapper::fileToTable(const char *filename) {
int res = 0; // DHT_from_file(dht_object, filename);
if (res != DHT_SUCCESS)
return res;
#ifdef DHT_STATISTICS
DAOSKV_print_statistics(daosKV_object);
#endif
return DHT_SUCCESS;
}
void DHT_Wrapper::printStatistics() {
int res;
res = DAOSKV_print_statistics(daosKV_object);
if (res != DHT_SUCCESS) {
// MPI ERROR ... WHAT TO DO NOW?
// RUNNING CIRCLES WHILE SCREAMING
}
}
LookupKey DHT_Wrapper::fuzzForDHT(int var_count, void *key, double dt) {
const auto c_zero_val = std::pow(10, AQUEOUS_EXP);
const Lookup_Keyelement dummy = {.0};
LookupKey vecFuzz(var_count + 1, dummy);
DHT_Rounder rounder;
int totals_i = 0;
// introduce fuzzing to allow more hits in DHT
// loop over every variable of grid cell
for (std::uint32_t i = 0; i < input_key_elements.size(); i++) {
if (input_key_elements[i] == DHT_KEY_INPUT_CUSTOM) {
continue;
}
double curr_key = ((double *)key)[input_key_elements[i]];
if (curr_key != 0) {
if (curr_key < c_zero_val &&
this->dht_prop_type_vector[i] == DHT_TYPE_DEFAULT) {
continue;
}
if (this->dht_prop_type_vector[i] == DHT_TYPE_TOTAL) {
curr_key -= base_totals[totals_i++];
}
vecFuzz[i] =
rounder.round(curr_key, dht_signif_vector[i],
this->dht_prop_type_vector[i] == DHT_TYPE_TOTAL);
}
}
// add timestep to the end of the key as double value
vecFuzz[var_count].fp_element = dt;
return vecFuzz;
}
void poet::DHT_Wrapper::SetSignifVector(std::vector<uint32_t> signif_vec) {
if (signif_vec.size() != this->key_count) {
throw std::runtime_error(
"Significant vector size mismatches count of key elements.");
}
this->dht_signif_vector = signif_vec;
}
} // namespace poet

View File

@ -0,0 +1,572 @@
/*
** Copyright (C) 2017-2021 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.
*/
#include "poet/DaosKeyValue.h"
#include <inttypes.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include <daos.h>
enum handleType
{
HANDLE_POOL,
HANDLE_CO,
};
static inline void
handle_share(DAOSKV *object, int type)
{
d_iov_t ghdl = {NULL, 0, 0};
int rc;
if (object->rank == 0)
{
/** fetch size of global handle */
if (type == HANDLE_POOL)
rc = daos_pool_local2global(object->poh, &ghdl);
else
rc = daos_cont_local2global(object->coh, &ghdl);
}
/** broadcast size of global handle to all peers */
MPI_Bcast(&ghdl.iov_buf_len, 1, MPI_UINT64_T, 0, object->communicator);
/** allocate buffer for global pool handle */
ghdl.iov_buf = malloc(ghdl.iov_buf_len);
ghdl.iov_len = ghdl.iov_buf_len;
if (object->rank == 0)
{
/** generate actual global handle to share with peer tasks */
if (type == HANDLE_POOL)
rc = daos_pool_local2global(object->poh, &ghdl);
else
rc = daos_cont_local2global(object->coh, &ghdl);
}
/** broadcast global handle to all peers */
MPI_Bcast(ghdl.iov_buf, ghdl.iov_len, MPI_BYTE, 0, object->communicator);
if (object->rank != 0)
{
/** unpack global handle */
if (type == HANDLE_POOL)
{
/* NB: Only pool_global2local are different */
rc = daos_pool_global2local(ghdl, &object->poh);
}
else
{
rc = daos_cont_global2local(object->poh, ghdl, &object->coh);
}
}
free(ghdl.iov_buf);
MPI_Barrier(object->communicator);
}
DAOSKV *DAOSKV_create(MPI_Comm comm)
{
DAOSKV *object;
int rc;
object = (DAOSKV *)malloc(sizeof(DAOSKV));
if (object == NULL)
return NULL;
// fill daos-object
object->communicator = comm;
object->read_misses = 0;
object->evictions = 0;
MPI_Comm_rank(comm, &object->rank);
MPI_Comm_size(comm, &object->comm_size);
// if set, initialize dht_stats
#ifdef DHT_STATISTICS
DAOSKV_stats *stats;
stats = (DAOSKV_stats *)malloc(sizeof(DAOSKV_stats));
if (stats == NULL)
return NULL;
object->stats = stats;
object->stats->writes_local = (int *)calloc(object->comm_size, sizeof(int));
object->stats->old_writes = 0;
object->stats->read_misses = 0;
object->stats->read_hits = 0;
object->stats->evictions = 0;
object->stats->w_access = 0;
object->stats->r_access = 0;
#endif
/** initialize the local DAOS stack */
if (daos_init() != 0)
return NULL;
/** Call connect on rank 0 only and broadcast handle to others */
if (object->rank == 0)
{
char *pool_name = getenv("DAOS_POOL");
if (pool_name == NULL)
printf("Pool label invalid \n");
if (daos_pool_connect(pool_name, NULL, DAOS_PC_RW, &object->poh,
NULL, NULL) != 0)
return NULL;
}
/** share pool handle with peer tasks */
handle_share(object, HANDLE_POOL);
/*
* Create and open container on rank 0 and share the handle.
*
* Alternatively, one could create the container outside of this program
* using the daos utility: daos cont create --pool=puuid
* and pass the uuid to the app.
*/
if (object->rank == 0)
{
/** create container */
if (getenv("DAOS_CONT") != NULL)
{
object->cont_label = getenv("DAOS_CONT");
}
else
{
object->cont_label = "Poet_Container";
}
/** check & open container if it already exist */
if (0 != daos_cont_open(object->poh, object->cont_label, DAOS_COO_RW, &object->coh, NULL, NULL))
{
/** create & open container*/
daos_cont_create_with_label(object->poh, object->cont_label, NULL, NULL, NULL);
/** open container */
if (daos_cont_open(object->poh, object->cont_label, DAOS_COO_RW, &object->coh, NULL, NULL) != 0)
return NULL;
}
}
/** share container handle with peer tasks */
handle_share(object, HANDLE_CO);
/** open object */
daos_obj_id_t oid;
oid.hi = 0;
oid.lo = 4;
daos_obj_generate_oid(object->coh, &oid, DAOS_OT_KV_HASHED, OC_SX, 0, 0);
daos_kv_open(object->coh, oid, DAOS_OO_RW, &object->oh, NULL);
return object;
}
int DAOSKV_free(DAOSKV *object)
{
MPI_Barrier(object->communicator);
if (daos_kv_close(object->oh, NULL) != 0)
return DAOS_ERROR;
if (daos_cont_close(object->coh, NULL) != 0)
return DAOS_ERROR;
if (object->rank == 0)
{
daos_cont_destroy(object->poh, object->cont_label, 0, NULL);
}
if (daos_pool_disconnect(object->poh, NULL) != 0)
return DAOS_ERROR;
if (daos_fini() != 0)
return DAOS_ERROR;
return DAOS_SUCCESS;
}
int DAOSKV_write(DAOSKV *object, void *key, int key_size, void *send_data, int send_size)
{
#ifdef DHT_STATISTICS
object->stats->w_access++;
#endif
// Turn key into a string
char *keyString[(key_size * 2) + 1];
keyToString(keyString, key, key_size);
int rc;
rc =daos_kv_put(object->oh, DAOS_TX_NONE, 0, keyString, send_size, send_data, NULL);
// No space left in storage
if (rc == -DER_NOSPACE && object->rank == 0)
{
trim_Space(object, 10, send_size, key_size);
}
if (rc != 0)
return DAOS_ERROR;
#ifdef DHT_STATISTICS
object->stats->writes_local[object->rank]++;
#endif
return DAOS_SUCCESS;
}
int DAOSKV_read(DAOSKV *object, void *key, int key_size, void *recv_data, int recv_size)
{
#ifdef DHT_STATISTICS
object->stats->r_access++;
#endif
// Turn key into a string
char *keyString[(key_size * 2) + 1];
keyToString(keyString, key, key_size);
daos_size_t size = recv_size;
int rc;
rc = daos_kv_get(object->oh, DAOS_TX_NONE, DAOS_COND_DKEY_FETCH, keyString, &size, recv_data, NULL);
if (rc == -DER_NONEXIST)
{
#ifdef DHT_STATISTICS
object->stats->read_misses += 1;
#endif
return DAOS_READ_MISS;
}
else if (rc != 0)
return DAOS_ERROR;
#ifdef DHT_STATISTICS
object->stats->read_hits += 1;
#endif
return DAOS_SUCCESS;
}
int DAOSKV_remove(DAOSKV *object, void *key, int key_size)
{
// Turn key into a string
char *keyString[(key_size * 2) + 1];
keyToString(keyString, key, key_size);
int rc;
if (daos_kv_remove(object->oh, DAOS_TX_NONE, 0, keyString, NULL) != 0)
return DAOS_ERROR;
return DAOS_SUCCESS;
}
int keyToString(char *output, void *key, int key_size)
{
int i;
int offset = 0;
for (i = 0; i < key_size; i++)
{
sprintf((char *)output + offset, "%02X", ((char *)key)[i]);
offset += 2;
}
output[offset++] = '\0';
return 0;
}
int enumerate_key(DAOSKV *object, int *total_nr, int key_size)
{
int actual_key_size = (key_size * 2) + 1;
char *buf;
daos_key_desc_t kds[5];
daos_anchor_t anchor = {0};
int key_nr = 0;
d_sg_list_t sgl;
d_iov_t sg_iov;
buf = malloc(actual_key_size);
d_iov_set(&sg_iov, buf, actual_key_size);
sgl.sg_nr = 1;
sgl.sg_nr_out = 0;
sgl.sg_iovs = &sg_iov;
while (!daos_anchor_is_eof(&anchor))
{
uint32_t nr = 5;
int rc;
memset(buf, 0, actual_key_size);
rc = daos_kv_list(object->oh, DAOS_TX_NONE, &nr, kds, &sgl, &anchor,
NULL);
//If there is no key, break the loop
if(buf[0] == '\0'){
break;
}
printf("Enumareted over dkey: %s\n", buf);
if (rc != 0)
{
printf("Error retrieving Key %d \n", rc);
return DAOS_ERROR;
}
if (nr == 0)
continue;
key_nr += nr;
}
*total_nr = key_nr;
return DAOS_SUCCESS;
}
int delete_n_entries(DAOSKV *object, int toDelete, int key_size)
{
int actual_key_size = (key_size * 2) + 1;
daos_handle_t th = DAOS_TX_NONE;
char *buf;
daos_key_desc_t kds[5];
daos_anchor_t anchor = {0};
d_sg_list_t sgl;
d_iov_t sg_iov;
int rc;
buf = malloc(actual_key_size);
d_iov_set(&sg_iov, buf, actual_key_size);
sgl.sg_nr = 1;
sgl.sg_nr_out = 0;
sgl.sg_iovs = &sg_iov;
memset(buf, 0, key_size);
/* allocate transaction */
rc = daos_tx_open(object->coh, &th, 0, NULL);
int key_nr = 0;
while (!daos_anchor_is_eof(&anchor) && key_nr < toDelete)
{
uint32_t nr = 5;
rc = daos_kv_list(object->oh, DAOS_TX_NONE, &nr, kds, &sgl, &anchor,
NULL);
//If there is no key, break the loop
if(buf[0] == '\0'){
break;
}
// Add delete of key to transaction th
printf("Delete dkey: %s\n", buf);
if (daos_kv_remove(object->oh, th, 0, buf, NULL) != 0)
printf("Delete n Key Error");
if (rc != 0)
{
printf("Error retrieving Key %d \n", rc);
return DAOS_ERROR;
}
if (nr == 0)
continue;
key_nr += nr;
}
// commit transaction, retry if failure
rc = daos_tx_commit(th, NULL);
if (rc)
{
printf("Commit error: %d\n", rc);
if (rc == -DER_TX_RESTART)
{
/* conflict with another transaction, try again */
rc = daos_tx_restart(th, NULL);
}
}
// free transaction resources
rc = daos_tx_close(th, NULL);
return DAOS_SUCCESS;
}
struct daos_space get_pool_size(DAOSKV *object)
{
int rc;
daos_pool_info_t pinfo = {0};
struct daos_pool_space *ps = &pinfo.pi_space;
// query only the space, replace with DPI_ALL for all infos
pinfo.pi_bits = DPI_SPACE;
rc = daos_pool_query(object->poh, NULL, &pinfo, NULL, NULL);
// size of storage
// printf("Total Size:%d\n", ps->ps_space.s_total[DAOS_MEDIA_SCM]+ps->ps_space.s_total[DAOS_MEDIA_NVME]);
// printf("Free Size:%d\n", ps->ps_space.s_free[DAOS_MEDIA_SCM]+ ps->ps_space.s_free[DAOS_MEDIA_NVME]);
return ps->ps_space;
}
int trim_Space(DAOSKV *object, float deletePercentage, int dataSize, int keySize)
{
// Get current usage of the storage space
struct daos_space space = get_pool_size(object);
long int total_size = space.s_total[DAOS_MEDIA_SCM] + space.s_total[DAOS_MEDIA_NVME];
// Estimate, total number of entries
int totalNumberOfEntries = total_size / (dataSize + keySize);
// Calculate how many keys to delete
int toDeleteEntries = totalNumberOfEntries * deletePercentage / 100;
delete_n_entries(object, toDeleteEntries, keySize);
return DAOS_SUCCESS;
}
int DAOSKV_print_statistics(DAOSKV *object)
{
#ifdef DHT_STATISTICS
int *written_buckets;
int *read_misses, sum_read_misses;
int *read_hits, sum_read_hits;
int *evictions, sum_evictions;
int sum_w_access, sum_r_access, *w_access, *r_access;
int rank;
MPI_Comm_rank(object->communicator, &rank);
// disable possible warning of unitialized variable, which is not the case
// since we're using MPI_Gather to obtain all values only on rank 0
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
// obtaining all values from all processes in the communicator
if (rank == 0)
read_misses = (int *)malloc(object->comm_size * sizeof(int));
if (MPI_Gather(&object->stats->read_misses, 1, MPI_INT, read_misses, 1,
MPI_INT, 0, object->communicator) != 0)
return DAOS_MPI_ERROR;
if (MPI_Reduce(&object->stats->read_misses, &sum_read_misses, 1, MPI_INT,
MPI_SUM, 0, object->communicator) != 0)
return DAOS_MPI_ERROR;
object->stats->read_misses = 0;
if (rank == 0)
read_hits = (int *)malloc(object->comm_size * sizeof(int));
if (MPI_Gather(&object->stats->read_hits, 1, MPI_INT, read_hits, 1,
MPI_INT, 0, object->communicator) != 0)
return DAOS_MPI_ERROR;
if (MPI_Reduce(&object->stats->read_hits, &sum_read_hits, 1, MPI_INT,
MPI_SUM, 0, object->communicator) != 0)
return DAOS_MPI_ERROR;
object->stats->read_hits = 0;
if (rank == 0)
evictions = (int *)malloc(object->comm_size * sizeof(int));
if (MPI_Gather(&object->stats->evictions, 1, MPI_INT, evictions, 1, MPI_INT, 0,
object->communicator) != 0)
return DAOS_MPI_ERROR;
if (MPI_Reduce(&object->stats->evictions, &sum_evictions, 1, MPI_INT, MPI_SUM,
0, object->communicator) != 0)
return DAOS_MPI_ERROR;
object->stats->evictions = 0;
if (rank == 0)
w_access = (int *)malloc(object->comm_size * sizeof(int));
if (MPI_Gather(&object->stats->w_access, 1, MPI_INT, w_access, 1, MPI_INT, 0,
object->communicator) != 0)
return DAOS_MPI_ERROR;
if (MPI_Reduce(&object->stats->w_access, &sum_w_access, 1, MPI_INT, MPI_SUM, 0,
object->communicator) != 0)
return DAOS_MPI_ERROR;
object->stats->w_access = 0;
if (rank == 0)
r_access = (int *)malloc(object->comm_size * sizeof(int));
if (MPI_Gather(&object->stats->r_access, 1, MPI_INT, r_access, 1, MPI_INT, 0,
object->communicator) != 0)
return DAOS_MPI_ERROR;
if (MPI_Reduce(&object->stats->r_access, &sum_r_access, 1, MPI_INT, MPI_SUM, 0,
object->communicator) != 0)
return DAOS_MPI_ERROR;
object->stats->r_access = 0;
if (rank == 0)
written_buckets = (int *)calloc(object->comm_size, sizeof(int));
if (MPI_Reduce(object->stats->writes_local, written_buckets, object->comm_size,
MPI_INT, MPI_SUM, 0, object->communicator) != 0)
return DAOS_MPI_ERROR;
if (rank == 0)
{ // only process with rank 0 will print out results as a
// object
int sum_written_buckets = 0;
for (int i = 0; i < object->comm_size; i++)
{
sum_written_buckets += written_buckets[i];
}
int members = 7;
int padsize = (members * 12) + 1;
char pad[padsize + 1];
memset(pad, '-', padsize * sizeof(char));
pad[padsize] = '\0';
printf("\n");
printf("%-35s||resets with every call of this function\n", " ");
printf("%-11s|%-11s|%-11s||%-11s|%-11s|%-11s|%-11s|%-11s\n", "rank", "occupied",
"free", "w_access", "r_access", "read misses", "read hits", "evictions");
printf("%s\n", pad);
for (int i = 0; i < object->comm_size; i++)
{
printf("%-11d|%-11d|%-11d||%-11d|%-11d|%-11d|%-11d|%-11d\n", i,
written_buckets[i], 0,
w_access[i], r_access[i], read_misses[i], read_hits[i], evictions[i]);
}
printf("%s\n", pad);
printf("%-11s|%-11d|%-11d||%-11d|%-11d|%-11d|%-11d|%-11d\n", "sum",
sum_written_buckets,
0,
sum_w_access, sum_r_access, sum_read_misses, sum_read_hits, sum_evictions);
printf("%s\n", pad);
printf("%s %d\n",
"new entries:", sum_written_buckets - object->stats->old_writes);
printf("\n");
fflush(stdout);
object->stats->old_writes = sum_written_buckets;
}
// enable warning again
#pragma GCC diagnostic pop
MPI_Barrier(object->communicator);
return DAOS_SUCCESS;
#endif
}

View File

@ -0,0 +1,597 @@
/*
** Copyright (C) 2017-2021 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.
*/
#include "poet/DaosKeyValue.h"
#include <inttypes.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include <daos.h>
enum handleType
{
HANDLE_POOL,
HANDLE_CO,
};
static inline void
handle_share(DAOSKV *object, int type)
{
d_iov_t ghdl = {NULL, 0, 0};
int rc;
if (object->rank == 0)
{
/** fetch size of global handle */
if (type == HANDLE_POOL)
rc = daos_pool_local2global(object->poh, &ghdl);
else
rc = daos_cont_local2global(object->coh, &ghdl);
}
/** broadcast size of global handle to all peers */
MPI_Bcast(&ghdl.iov_buf_len, 1, MPI_UINT64_T, 0, object->communicator);
/** allocate buffer for global pool handle */
ghdl.iov_buf = malloc(ghdl.iov_buf_len);
ghdl.iov_len = ghdl.iov_buf_len;
if (object->rank == 0)
{
/** generate actual global handle to share with peer tasks */
if (type == HANDLE_POOL)
rc = daos_pool_local2global(object->poh, &ghdl);
else
rc = daos_cont_local2global(object->coh, &ghdl);
}
/** broadcast global handle to all peers */
MPI_Bcast(ghdl.iov_buf, ghdl.iov_len, MPI_BYTE, 0, object->communicator);
if (object->rank != 0)
{
/** unpack global handle */
if (type == HANDLE_POOL)
{
/* NB: Only pool_global2local are different */
rc = daos_pool_global2local(ghdl, &object->poh);
}
else
{
rc = daos_cont_global2local(object->poh, ghdl, &object->coh);
}
}
free(ghdl.iov_buf);
MPI_Barrier(object->communicator);
}
DAOSKV *DAOSKV_create(MPI_Comm comm)
{
DAOSKV *object;
int rc;
object = (DAOSKV *)malloc(sizeof(DAOSKV));
if (object == NULL)
return NULL;
// fill daos-object
object->communicator = comm;
object->read_misses = 0;
object->evictions = 0;
MPI_Comm_rank(comm, &object->rank);
MPI_Comm_size(comm, &object->comm_size);
// if set, initialize dht_stats
#ifdef DHT_STATISTICS
DAOSKV_stats *stats;
stats = (DAOSKV_stats *)malloc(sizeof(DAOSKV_stats));
if (stats == NULL)
return NULL;
object->stats = stats;
object->stats->writes_local = (int *)calloc(object->comm_size, sizeof(int));
object->stats->old_writes = 0;
object->stats->read_misses = 0;
object->stats->read_hits = 0;
object->stats->evictions = 0;
object->stats->w_access = 0;
object->stats->r_access = 0;
#endif
/** initialize the local DAOS stack */
if (daos_init() != 0)
return NULL;
/** Call connect on rank 0 only and broadcast handle to others */
if (object->rank == 0)
{
char *pool_name = getenv("DAOS_POOL");
if (pool_name == NULL)
printf("Pool label invalid \n");
if (daos_pool_connect(pool_name, NULL, DAOS_PC_RW, &object->poh,
NULL, NULL) != 0)
return NULL;
}
/** share pool handle with peer tasks */
handle_share(object, HANDLE_POOL);
/*
* Create and open container on rank 0 and share the handle.
*
* Alternatively, one could create the container outside of this program
* using the daos utility: daos cont create --pool=puuid
* and pass the uuid to the app.
*/
if (object->rank == 0)
{
/** create container */
if (getenv("DAOS_CONT") != NULL)
{
object->cont_label = getenv("DAOS_CONT");
}
else
{
object->cont_label = "Poet_Container";
}
/** check & open container if it already exist */
if (0 != daos_cont_open(object->poh, object->cont_label, DAOS_COO_RW, &object->coh, NULL, NULL))
{
/** create & open container*/
daos_cont_create_with_label(object->poh, object->cont_label, NULL, NULL, NULL);
/** open container */
if (daos_cont_open(object->poh, object->cont_label, DAOS_COO_RW, &object->coh, NULL, NULL) != 0)
return NULL;
}
}
/** share container handle with peer tasks */
handle_share(object, HANDLE_CO);
/** open object */
daos_obj_id_t oid;
oid.hi = 0;
oid.lo = 4;
daos_obj_generate_oid(object->coh, &oid, 0, OC_SX, 0, 0);
daos_obj_open(object->coh, oid, DAOS_OO_RW, &object->oh, NULL);
return object;
}
int DAOSKV_free(DAOSKV *object)
{
MPI_Barrier(object->communicator);
if (daos_obj_close(object->oh, NULL) != 0)
return DAOS_ERROR;
if (daos_cont_close(object->coh, NULL) != 0)
return DAOS_ERROR;
if (object->rank == 0)
{
daos_cont_destroy(object->poh, object->cont_label, 0, NULL);
}
if (daos_pool_disconnect(object->poh, NULL) != 0)
return DAOS_ERROR;
if (daos_fini() != 0)
return DAOS_ERROR;
return DAOS_SUCCESS;
}
int DAOSKV_write(DAOSKV *object, void *key, int key_size, void *send_data, int send_size)
{
int rc;
#ifdef DHT_STATISTICS
object->stats->w_access++;
#endif
d_iov_t dkey;
d_sg_list_t sgl;
d_iov_t sg_iov;
daos_iod_t iod;
// set dkey
d_iov_set(&dkey, key, key_size);
d_iov_set(&sg_iov, send_data, send_size);
sgl.sg_nr = 1;
sgl.sg_nr_out = 0;
sgl.sg_iovs = &sg_iov;
// set akey
// Here for all the same value since dkey differ
// Maybe other way around is better? Or a mix?
int akey = 1;
d_iov_set(&iod.iod_name, &akey, sizeof(int));
iod.iod_nr = 1; /** has to be 1 for single value */
iod.iod_size = send_size; /** size of the single value */
iod.iod_recxs = NULL; /** recx is ignored for single value */
iod.iod_type = DAOS_IOD_SINGLE; /** value type of the akey */
rc = daos_obj_update(object->oh, DAOS_TX_NONE, 0, &dkey, 1, &iod, &sgl,
NULL);
// No space left in storage
if (rc == -DER_NOSPACE && object->rank == 0)
{
trim_Space(object, 10, send_size, key_size);
}
if (rc != 0)
return DAOS_ERROR;
#ifdef DHT_STATISTICS
object->stats->writes_local[object->rank]++;
#endif
return DAOS_SUCCESS;
}
int DAOSKV_read(DAOSKV *object, void *key, int key_size, void *recv_data, int recv_size)
{
int rc;
#ifdef DHT_STATISTICS
object->stats->r_access++;
#endif
d_iov_t dkey;
d_sg_list_t sgl;
d_iov_t sg_iov;
daos_iod_t iod;
// set dkey
d_iov_set(&dkey, key, key_size);
d_iov_set(&sg_iov, recv_data, recv_size);
sgl.sg_nr = 1;
sgl.sg_nr_out = 0;
sgl.sg_iovs = &sg_iov;
// set akey
// Here for all the same value since dkey differ
// Maybe other way around is better? Or a mix?
int akey = 1;
d_iov_set(&iod.iod_name, &akey, sizeof(int));
iod.iod_nr = 1; /** 1 for single value */
iod.iod_size = DAOS_REC_ANY; /** size of the single value, set to DAOS_REC_ANY to check if key was written or not */
iod.iod_recxs = NULL; /** recx is ignored for single value */
iod.iod_type = DAOS_IOD_SINGLE; /** value type of the akey */
/** fetch a dkey */
rc = daos_obj_fetch(object->oh, DAOS_TX_NONE, 0, &dkey, 1, &iod, &sgl,
NULL, NULL);
if (rc != 0)
{
return DAOS_ERROR;
}
if (iod.iod_size == 0)
{
#ifdef DHT_STATISTICS
object->stats->read_misses += 1;
#endif
return DAOS_READ_MISS;
}
#ifdef DHT_STATISTICS
object->stats->read_hits += 1;
#endif
return DAOS_SUCCESS;
}
int DAOSKV_remove(DAOSKV *object, void *key, int key_size)
{
d_iov_t dkey;
// set dkey
d_iov_set(&dkey, key, key_size);
if (daos_obj_punch_dkeys(object->oh, DAOS_TX_NONE, 0, 1, &dkey, NULL) != 0)
return DAOS_ERROR;
return DAOS_SUCCESS;
}
int enumerate_key(DAOSKV *object, int *total_nr, int key_size)
{
char *buf;
daos_key_desc_t kds[5];
daos_anchor_t anchor = {0};
d_sg_list_t sgl;
d_iov_t sg_iov;
int key_nr = 0;
int rc;
buf = malloc(key_size);
d_iov_set(&sg_iov, buf, key_size);
sgl.sg_nr = 1;
sgl.sg_nr_out = 0;
sgl.sg_iovs = &sg_iov;
while (!daos_anchor_is_eof(&anchor))
{
uint32_t nr = 5;
memset(buf, 0, key_size);
rc = daos_obj_list_dkey(object->oh, DAOS_TX_NONE, &nr, kds,
&sgl, &anchor, NULL);
//If there is no key, break the loop
if(buf[0] == '\0'){
break;
}
printf("Enumareted over dkey: %s\n", buf);
if (rc != 0)
{
printf("Error retrieving Key %d \n", rc);
return DAOS_ERROR;
}
if (nr == 0)
continue;
key_nr += nr;
}
*total_nr = key_nr;
return DAOS_SUCCESS;
}
int delete_n_entries(DAOSKV *object, int toDelete, int key_size)
{
daos_handle_t th = DAOS_TX_NONE;
char *buf;
daos_key_desc_t kds[5];
daos_anchor_t anchor = {0};
d_sg_list_t sgl;
d_iov_t sg_iov;
int rc;
buf = malloc(key_size);
d_iov_set(&sg_iov, buf, key_size);
sgl.sg_nr = 1;
sgl.sg_nr_out = 0;
sgl.sg_iovs = &sg_iov;
memset(buf, 0, key_size);
/* allocate transaction */
rc = daos_tx_open(object->coh, &th, 0, NULL);
int key_nr = 0;
while (!daos_anchor_is_eof(&anchor) && key_nr < toDelete)
{
uint32_t nr = 5;
rc = daos_obj_list_dkey(object->oh, DAOS_TX_NONE, &nr, kds,
&sgl, &anchor, NULL);
//If there is no key, break the loop
if(buf[0] == '\0'){
break;
}
// Add delete of key to transaction th
printf("Delete dkey: %s\n", buf);
d_iov_t dkey;
// set dkey
d_iov_set(&dkey, buf, key_size);
if (daos_obj_punch_dkeys(object->oh, th, 0, 1, &dkey, NULL) != 0)
printf("Delete n Key Error");
if (rc != 0)
{
printf("Error retrieving Key %d \n", rc);
return DAOS_ERROR;
}
if (nr == 0)
continue;
key_nr += nr;
}
// commit transaction, retry if failure
rc = daos_tx_commit(th, NULL);
if (rc)
{
printf("Commit error: %d\n", rc);
if (rc == -DER_TX_RESTART)
{
/* conflict with another transaction, try again */
rc = daos_tx_restart(th, NULL);
}
}
// free transaction resources
rc = daos_tx_close(th, NULL);
return DAOS_SUCCESS;
}
struct daos_space get_pool_size(DAOSKV *object)
{
int rc;
daos_pool_info_t pinfo = {0};
struct daos_pool_space *ps = &pinfo.pi_space;
// query only the space, replace with DPI_ALL for all infos
pinfo.pi_bits = DPI_SPACE;
rc = daos_pool_query(object->poh, NULL, &pinfo, NULL, NULL);
// size of storage
// printf("Total Size:%d\n", ps->ps_space.s_total[DAOS_MEDIA_SCM]+ps->ps_space.s_total[DAOS_MEDIA_NVME]);
// printf("Free Size:%d\n", ps->ps_space.s_free[DAOS_MEDIA_SCM]+ ps->ps_space.s_free[DAOS_MEDIA_NVME]);
return ps->ps_space;
}
int trim_Space(DAOSKV *object, float deletePercentage, int dataSize, int keySize)
{
// Get current usage of the storage space
struct daos_space space = get_pool_size(object);
long int total_size = space.s_total[DAOS_MEDIA_SCM] + space.s_total[DAOS_MEDIA_NVME];
// Estimate, total number of entries
int totalNumberOfEntries = total_size / (dataSize + keySize);
// Calculate how many keys to delete
int toDeleteEntries = totalNumberOfEntries * deletePercentage / 100;
delete_n_entries(object, toDeleteEntries, keySize);
return DAOS_SUCCESS;
}
int DAOSKV_print_statistics(DAOSKV *object)
{
#ifdef DHT_STATISTICS
int *written_buckets;
int *read_misses, sum_read_misses;
int *read_hits, sum_read_hits;
int *evictions, sum_evictions;
int sum_w_access, sum_r_access, *w_access, *r_access;
int rank;
MPI_Comm_rank(object->communicator, &rank);
// disable possible warning of unitialized variable, which is not the case
// since we're using MPI_Gather to obtain all values only on rank 0
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
// obtaining all values from all processes in the communicator
if (rank == 0)
read_misses = (int *)malloc(object->comm_size * sizeof(int));
if (MPI_Gather(&object->stats->read_misses, 1, MPI_INT, read_misses, 1,
MPI_INT, 0, object->communicator) != 0)
return DAOS_MPI_ERROR;
if (MPI_Reduce(&object->stats->read_misses, &sum_read_misses, 1, MPI_INT,
MPI_SUM, 0, object->communicator) != 0)
return DAOS_MPI_ERROR;
object->stats->read_misses = 0;
if (rank == 0)
read_hits = (int *)malloc(object->comm_size * sizeof(int));
if (MPI_Gather(&object->stats->read_hits, 1, MPI_INT, read_hits, 1,
MPI_INT, 0, object->communicator) != 0)
return DAOS_MPI_ERROR;
if (MPI_Reduce(&object->stats->read_hits, &sum_read_hits, 1, MPI_INT,
MPI_SUM, 0, object->communicator) != 0)
return DAOS_MPI_ERROR;
object->stats->read_hits = 0;
if (rank == 0)
evictions = (int *)malloc(object->comm_size * sizeof(int));
if (MPI_Gather(&object->stats->evictions, 1, MPI_INT, evictions, 1, MPI_INT, 0,
object->communicator) != 0)
return DAOS_MPI_ERROR;
if (MPI_Reduce(&object->stats->evictions, &sum_evictions, 1, MPI_INT, MPI_SUM,
0, object->communicator) != 0)
return DAOS_MPI_ERROR;
object->stats->evictions = 0;
if (rank == 0)
w_access = (int *)malloc(object->comm_size * sizeof(int));
if (MPI_Gather(&object->stats->w_access, 1, MPI_INT, w_access, 1, MPI_INT, 0,
object->communicator) != 0)
return DAOS_MPI_ERROR;
if (MPI_Reduce(&object->stats->w_access, &sum_w_access, 1, MPI_INT, MPI_SUM, 0,
object->communicator) != 0)
return DAOS_MPI_ERROR;
object->stats->w_access = 0;
if (rank == 0)
r_access = (int *)malloc(object->comm_size * sizeof(int));
if (MPI_Gather(&object->stats->r_access, 1, MPI_INT, r_access, 1, MPI_INT, 0,
object->communicator) != 0)
return DAOS_MPI_ERROR;
if (MPI_Reduce(&object->stats->r_access, &sum_r_access, 1, MPI_INT, MPI_SUM, 0,
object->communicator) != 0)
return DAOS_MPI_ERROR;
object->stats->r_access = 0;
if (rank == 0)
written_buckets = (int *)calloc(object->comm_size, sizeof(int));
if (MPI_Reduce(object->stats->writes_local, written_buckets, object->comm_size,
MPI_INT, MPI_SUM, 0, object->communicator) != 0)
return DAOS_MPI_ERROR;
if (rank == 0)
{ // only process with rank 0 will print out results as a
// object
int sum_written_buckets = 0;
for (int i = 0; i < object->comm_size; i++)
{
sum_written_buckets += written_buckets[i];
}
int members = 7;
int padsize = (members * 12) + 1;
char pad[padsize + 1];
memset(pad, '-', padsize * sizeof(char));
pad[padsize] = '\0';
printf("\n");
printf("%-35s||resets with every call of this function\n", " ");
printf("%-11s|%-11s|%-11s||%-11s|%-11s|%-11s|%-11s|%-11s\n", "rank", "occupied",
"free", "w_access", "r_access", "read misses", "read hits", "evictions");
printf("%s\n", pad);
for (int i = 0; i < object->comm_size; i++)
{
printf("%-11d|%-11d|%-11d||%-11d|%-11d|%-11d|%-11d|%-11d\n", i,
written_buckets[i], 0,
w_access[i], r_access[i], read_misses[i], read_hits[i], evictions[i]);
}
printf("%s\n", pad);
printf("%-11s|%-11d|%-11d||%-11d|%-11d|%-11d|%-11d|%-11d\n", "sum",
sum_written_buckets,
0,
sum_w_access, sum_r_access, sum_read_misses, sum_read_hits, sum_evictions);
printf("%s\n", pad);
printf("%s %d\n",
"new entries:", sum_written_buckets - object->stats->old_writes);
printf("\n");
fflush(stdout);
object->stats->old_writes = sum_written_buckets;
}
// enable warning again
#pragma GCC diagnostic pop
MPI_Barrier(object->communicator);
return DAOS_SUCCESS;
#endif
}

View File

@ -1,215 +0,0 @@
// Time-stamp: "Last modified 2023-07-21 17:22:00 mluebke"
/*
** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
** Potsdam)
**
** Copyright (C) 2018-2021 Marco De Lucia (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 "poet/DHT_Wrapper.hpp"
#include "poet/DHT_Types.hpp"
#include "poet/HashFunctions.hpp"
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <cstring>
#include <iostream>
#include <math.h>
#include <stdexcept>
#include <vector>
using namespace poet;
using namespace std;
inline DHT_SCNotation round_key_element(double value, std::uint32_t signif) {
std::int8_t exp =
static_cast<std::int8_t>(std::floor(std::log10(std::fabs(value))));
std::int64_t significant =
static_cast<std::int64_t>(value * std::pow(10, signif - exp - 1));
DHT_SCNotation elem;
elem.exp = exp;
elem.significant = significant;
return elem;
}
DHT_Wrapper::DHT_Wrapper(MPI_Comm dht_comm, uint32_t dht_size,
const std::vector<std::uint32_t> &key_indices,
uint32_t data_count)
: key_count(key_indices.size()), data_count(data_count),
input_key_elements(key_indices) {
// initialize DHT object
uint32_t key_size = (key_count + 1) * 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::Murmur2_64A);
this->dht_signif_vector.resize(key_size, DHT_KEY_SIGNIF_DEFAULT);
this->dht_prop_type_vector.resize(key_count, DHT_TYPE_DEFAULT);
this->dht_prop_type_vector[0] = DHT_TYPE_TOTAL;
this->dht_prop_type_vector[1] = DHT_TYPE_TOTAL;
this->dht_prop_type_vector[2] = DHT_TYPE_CHARGE;
}
DHT_Wrapper::~DHT_Wrapper() {
// free DHT
DHT_free(dht_object, NULL, NULL);
}
auto DHT_Wrapper::checkDHT(int length, double dt,
const std::vector<double> &work_package,
std::vector<std::uint32_t> &curr_mapping)
-> const poet::DHT_ResultObject & {
dht_results.length = length;
dht_results.keys.resize(length);
dht_results.results.resize(length);
dht_results.needPhreeqc.resize(length);
std::vector<std::uint32_t> new_mapping;
// loop over every grid cell contained in work package
for (int i = 0; i < length; i++) {
// point to current grid cell
void *key = (void *)&(work_package[i * this->data_count]);
auto &data = dht_results.results[i];
auto &key_vector = dht_results.keys[i];
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
int res = DHT_read(this->dht_object, key_vector.data(), data.data());
switch (res) {
case DHT_SUCCESS:
dht_results.needPhreeqc[i] = false;
this->dht_hits++;
break;
case DHT_READ_MISS:
dht_results.needPhreeqc[i] = true;
new_mapping.push_back(curr_mapping[i]);
this->dht_miss++;
break;
}
}
curr_mapping = std::move(new_mapping);
return dht_results;
}
void DHT_Wrapper::fillDHT(int length, const std::vector<double> &work_package) {
// loop over every grid cell contained in work package
for (int i = 0; i < length; i++) {
// If true grid cell was simulated, needs to be inserted into dht
if (dht_results.needPhreeqc[i]) {
const auto &key = dht_results.keys[i];
void *data = (void *)&(work_package[i * this->data_count]);
// fuzz data (round, logarithm etc.)
// insert simulated data with fuzzed key into DHT
int res = DHT_write(this->dht_object, (void *)key.data(), data);
// if data was successfully written ...
if ((res != DHT_SUCCESS) && (res == DHT_WRITE_SUCCESS_WITH_EVICTION)) {
dht_evictions++;
}
}
}
}
void DHT_Wrapper::resultsToWP(std::vector<double> &work_package) {
for (int i = 0; i < dht_results.length; i++) {
if (!dht_results.needPhreeqc[i]) {
std::copy(dht_results.results[i].begin(), dht_results.results[i].end(),
work_package.begin() + (data_count * i));
}
}
}
int DHT_Wrapper::tableToFile(const char *filename) {
int res = DHT_to_file(dht_object, filename);
return res;
}
int DHT_Wrapper::fileToTable(const char *filename) {
int res = DHT_from_file(dht_object, filename);
if (res != DHT_SUCCESS)
return res;
#ifdef DHT_STATISTICS
DHT_print_statistics(dht_object);
#endif
return DHT_SUCCESS;
}
void DHT_Wrapper::printStatistics() {
int res;
res = DHT_print_statistics(dht_object);
if (res != DHT_SUCCESS) {
// MPI ERROR ... WHAT TO DO NOW?
// RUNNING CIRCLES WHILE SCREAMING
}
}
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 + 1);
std::memset(&vecFuzz[0], 0, sizeof(DHT_Keyelement) * var_count);
int totals_i = 0;
// introduce fuzzing to allow more hits in DHT
// loop over every variable of grid cell
for (std::uint32_t i = 0; i < input_key_elements.size(); i++) {
double curr_key = ((double *)key)[input_key_elements[i]];
if (curr_key != 0) {
if (curr_key < zero_val &&
this->dht_prop_type_vector[i] == DHT_TYPE_DEFAULT) {
continue;
}
if (this->dht_prop_type_vector[i] == DHT_TYPE_TOTAL) {
curr_key -= base_totals[totals_i++];
}
vecFuzz[i].sc_notation =
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
vecFuzz[var_count].fp_elemet = dt;
return vecFuzz;
}
void poet::DHT_Wrapper::SetSignifVector(std::vector<uint32_t> signif_vec) {
if (signif_vec.size() != this->key_count) {
throw std::runtime_error(
"Significant vector size mismatches count of key elements.");
}
this->dht_signif_vector = signif_vec;
}

153
src/ChemistryModule/MasterFunctions.cpp Normal file → Executable file
View File

@ -2,6 +2,7 @@
#include "poet/ChemistryModule.hpp"
#include <algorithm>
#include <cstddef>
#include <cstdint>
#include <mpi.h>
#include <stdexcept>
@ -62,19 +63,136 @@ std::vector<double> poet::ChemistryModule::GetWorkerIdleTimings() const {
std::vector<uint32_t> poet::ChemistryModule::GetWorkerDHTHits() const {
int type = CHEM_PERF;
MPI_Bcast(&type, 1, MPI_INT, 0, this->group_comm);
return MasterGatherWorkerMetrics(WORKER_DHT_HITS);
}
std::vector<uint32_t> poet::ChemistryModule::GetWorkerDHTMiss() const {
int type = CHEM_PERF;
type = WORKER_DHT_HITS;
MPI_Bcast(&type, 1, MPI_INT, 0, this->group_comm);
return MasterGatherWorkerMetrics(WORKER_DHT_MISS);
MPI_Status probe;
MPI_Probe(MPI_ANY_SOURCE, WORKER_DHT_HITS, this->group_comm, &probe);
int count;
MPI_Get_count(&probe, MPI_UINT32_T, &count);
std::vector<uint32_t> ret(count);
MPI_Recv(ret.data(), count, MPI_UINT32_T, probe.MPI_SOURCE, WORKER_DHT_HITS,
this->group_comm, NULL);
return ret;
}
std::vector<uint32_t> poet::ChemistryModule::GetWorkerDHTEvictions() const {
int type = CHEM_PERF;
MPI_Bcast(&type, 1, MPI_INT, 0, this->group_comm);
return MasterGatherWorkerMetrics(WORKER_DHT_EVICTIONS);
type = WORKER_DHT_EVICTIONS;
MPI_Bcast(&type, 1, MPI_INT, 0, this->group_comm);
MPI_Status probe;
MPI_Probe(MPI_ANY_SOURCE, WORKER_DHT_EVICTIONS, this->group_comm, &probe);
int count;
MPI_Get_count(&probe, MPI_UINT32_T, &count);
std::vector<uint32_t> ret(count);
MPI_Recv(ret.data(), count, MPI_UINT32_T, probe.MPI_SOURCE, WORKER_DHT_EVICTIONS,
this->group_comm, NULL);
return ret;
}
std::vector<double>
poet::ChemistryModule::GetWorkerInterpolationWriteTimings() const {
int type = CHEM_PERF;
MPI_Bcast(&type, 1, MPI_INT, 0, this->group_comm);
return MasterGatherWorkerTimings(WORKER_IP_WRITE);
}
std::vector<double>
poet::ChemistryModule::GetWorkerInterpolationReadTimings() const {
int type = CHEM_PERF;
MPI_Bcast(&type, 1, MPI_INT, 0, this->group_comm);
return MasterGatherWorkerTimings(WORKER_IP_READ);
}
std::vector<double>
poet::ChemistryModule::GetWorkerInterpolationGatherTimings() const {
int type = CHEM_PERF;
MPI_Bcast(&type, 1, MPI_INT, 0, this->group_comm);
return MasterGatherWorkerTimings(WORKER_IP_GATHER);
}
std::vector<double>
poet::ChemistryModule::GetWorkerInterpolationFunctionCallTimings() const {
int type = CHEM_PERF;
MPI_Bcast(&type, 1, MPI_INT, 0, this->group_comm);
return MasterGatherWorkerTimings(WORKER_IP_FC);
}
std::vector<uint32_t>
poet::ChemistryModule::GetWorkerInterpolationCalls() const {
int type = CHEM_PERF;
MPI_Bcast(&type, 1, MPI_INT, 0, this->group_comm);
type = WORKER_IP_CALLS;
MPI_Bcast(&type, 1, MPI_INT, 0, this->group_comm);
MPI_Status probe;
MPI_Probe(MPI_ANY_SOURCE, WORKER_IP_CALLS, this->group_comm, &probe);
int count;
MPI_Get_count(&probe, MPI_UINT32_T, &count);
std::vector<uint32_t> ret(count);
MPI_Recv(ret.data(), count, MPI_UINT32_T, probe.MPI_SOURCE, WORKER_IP_CALLS,
this->group_comm, NULL);
return ret;
}
std::vector<uint32_t> poet::ChemistryModule::GetWorkerPHTCacheHits() const {
int type = CHEM_PERF;
MPI_Bcast(&type, 1, MPI_INT, 0, this->group_comm);
type = WORKER_PHT_CACHE_HITS;
MPI_Bcast(&type, 1, MPI_INT, 0, this->group_comm);
MPI_Status probe;
MPI_Probe(MPI_ANY_SOURCE, type, this->group_comm, &probe);
int count;
MPI_Get_count(&probe, MPI_UINT32_T, &count);
std::vector<uint32_t> ret(count);
MPI_Recv(ret.data(), count, MPI_UINT32_T, probe.MPI_SOURCE, type,
this->group_comm, NULL);
return ret;
}
inline std::vector<double> shuffleField(const std::vector<double> &in_field,
uint32_t size_per_prop,
uint32_t prop_count,
uint32_t wp_count) {
std::vector<double> out_buffer(in_field.size());
uint32_t write_i = 0;
for (uint32_t i = 0; i < wp_count; i++) {
for (uint32_t j = i; j < size_per_prop; j += wp_count) {
for (uint32_t k = 0; k < prop_count; k++) {
out_buffer[(write_i * prop_count) + k] =
in_field[(k * size_per_prop) + j];
}
write_i++;
}
}
return out_buffer;
}
inline void unshuffleField(const std::vector<double> &in_buffer,
uint32_t size_per_prop, uint32_t prop_count,
uint32_t wp_count, std::vector<double> &out_field) {
uint32_t read_i = 0;
for (uint32_t i = 0; i < wp_count; i++) {
for (uint32_t j = i; j < size_per_prop; j += wp_count) {
for (uint32_t k = 0; k < prop_count; k++) {
out_field[(k * size_per_prop) + j] =
in_buffer[(read_i * prop_count) + k];
}
read_i++;
}
}
}
inline void printProgressbar(int count_pkgs, int n_wp, int barWidth = 70) {
@ -190,12 +308,15 @@ inline void poet::ChemistryModule::MasterRecvPkgs(worker_list_t &w_list,
}
void poet::ChemistryModule::RunCells() {
double start_t{MPI_Wtime()};
if (this->is_sequential) {
MasterRunSequential();
return;
}
MasterRunParallel();
double end_t{MPI_Wtime()};
this->chem_t += end_t - start_t;
}
void poet::ChemistryModule::MasterRunSequential() {
@ -211,7 +332,6 @@ void poet::ChemistryModule::MasterRunSequential() {
void poet::ChemistryModule::MasterRunParallel() {
/* declare most of the needed variables here */
double chem_a, chem_b;
double seq_a, seq_b, seq_c, seq_d;
double worker_chemistry_a, worker_chemistry_b;
double sim_e_chemistry, sim_f_chemistry;
@ -227,9 +347,6 @@ void poet::ChemistryModule::MasterRunParallel() {
double dt = this->PhreeqcRM::GetTimeStep();
static uint32_t iteration = 0;
/* start time measurement of whole chemistry simulation */
chem_a = MPI_Wtime();
/* start time measurement of sequential part */
seq_a = MPI_Wtime();
@ -308,8 +425,6 @@ void poet::ChemistryModule::MasterRunParallel() {
for (int i = 1; i < this->comm_size; i++) {
MPI_Send(NULL, 0, MPI_DOUBLE, i, LOOP_END, this->group_comm);
}
chem_b = MPI_Wtime();
this->chem_t += chem_b - chem_a;
this->simtime += dt;
iteration++;
@ -324,7 +439,8 @@ std::vector<uint32_t>
poet::ChemistryModule::CalculateWPSizesVector(uint32_t n_cells,
uint32_t wp_size) const {
bool mod_pkgs = (n_cells % wp_size) != 0;
uint32_t n_packages = (uint32_t)(n_cells / wp_size) + mod_pkgs;
uint32_t n_packages =
(uint32_t)(n_cells / wp_size) + static_cast<int>(mod_pkgs);
std::vector<uint32_t> wp_sizes_vector(n_packages, 0);
@ -334,12 +450,3 @@ poet::ChemistryModule::CalculateWPSizesVector(uint32_t n_cells,
return wp_sizes_vector;
}
void poet::ChemistryModule::setProgressBarPrintout(bool enabled) {
if (is_master) {
int type = CHEM_PROGRESSBAR;
ChemBCast(&type, 1, MPI_INT);
ChemBCast(&enabled, 1, MPI_CXX_BOOL);
}
this->print_progessbar = enabled;
}

View File

@ -0,0 +1,21 @@
file(GLOB surrogate_models_SRC
CONFIGURE_DEPENDS
"*.cpp" "*.c")
find_library(MATH_LIBRARY m)
add_library(surrogate_models ${surrogate_models_SRC})
target_include_directories(surrogate_models PUBLIC ${PROJECT_SOURCE_DIR}/include)
target_link_libraries(surrogate_models PUBLIC
MPI::MPI_CXX ${MATH_LIBRARY})
target_compile_definitions(surrogate_models PUBLIC OMPI_SKIP_MPICXX)
option(POET_PHT_ADDITIONAL_INFO "Enables additional information in the PHT" OFF)
if (POET_PHT_ADDITIONAL_INFO)
target_compile_definitions(surrogate_models PRIVATE POET_PHT_ADD)
endif()
if(POET_DHT_DEBUG)
target_compile_definitions(surrogate_models PRIVATE DHT_STATISTICS)
endif()

View File

@ -1,3 +1,4 @@
/// Time-stamp: "Last modified 2023-06-28 15:58:19 mluebke"
/*
** Copyright (C) 2017-2021 Max Luebke (University of Potsdam)
**
@ -15,10 +16,12 @@
** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
#include <mpi.h>
#include <poet/DHT.h>
#include <inttypes.h>
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
@ -52,7 +55,8 @@ static int read_flag(char flag_byte) {
}
DHT *DHT_create(MPI_Comm comm, uint64_t size, unsigned int data_size,
unsigned int key_size, uint64_t (*hash_func)(int, const void *)) {
unsigned int key_size,
uint64_t (*hash_func)(int, const void *)) {
DHT *object;
MPI_Win window;
void *mem_alloc;
@ -61,17 +65,20 @@ DHT *DHT_create(MPI_Comm comm, uint64_t size, unsigned int data_size,
// calculate how much bytes for the index are needed to address count of
// buckets per process
index_bytes = (int)ceil(log2(size));
if (index_bytes % 8 != 0) index_bytes = index_bytes + (8 - (index_bytes % 8));
if (index_bytes % 8 != 0)
index_bytes = index_bytes + (8 - (index_bytes % 8));
// allocate memory for dht-object
object = (DHT *)malloc(sizeof(DHT));
if (object == NULL) return NULL;
if (object == NULL)
return NULL;
// every memory allocation has 1 additional byte for flags etc.
if (MPI_Alloc_mem(size * (1 + data_size + key_size), MPI_INFO_NULL,
&mem_alloc) != 0)
return NULL;
if (MPI_Comm_size(comm, &comm_size) != 0) return NULL;
if (MPI_Comm_size(comm, &comm_size) != 0)
return NULL;
// since MPI_Alloc_mem doesn't provide memory allocation with the memory set
// to zero, we're doing this here
@ -104,7 +111,8 @@ DHT *DHT_create(MPI_Comm comm, uint64_t size, unsigned int data_size,
DHT_stats *stats;
stats = (DHT_stats *)malloc(sizeof(DHT_stats));
if (stats == NULL) return NULL;
if (stats == NULL)
return NULL;
object->stats = stats;
object->stats->writes_local = (int *)calloc(comm_size, sizeof(int));
@ -118,7 +126,106 @@ DHT *DHT_create(MPI_Comm comm, uint64_t size, unsigned int data_size,
return object;
}
int DHT_write(DHT *table, void *send_key, void *send_data) {
void DHT_set_accumulate_callback(DHT *table,
int (*callback_func)(int, void *, int,
void *)) {
table->accumulate_callback = callback_func;
}
int DHT_write_accumulate(DHT *table, const void *send_key, int data_size,
void *send_data, uint32_t *proc, uint32_t *index,
int *callback_ret) {
unsigned int dest_rank, i;
int result = DHT_SUCCESS;
#ifdef DHT_STATISTICS
table->stats->w_access++;
#endif
// determine destination rank and index by hash of key
determine_dest(table->hash_func(table->key_size, send_key), table->comm_size,
table->table_size, &dest_rank, table->index,
table->index_count);
// concatenating key with data to write entry to DHT
set_flag((char *)table->send_entry);
memcpy((char *)table->send_entry + 1, (char *)send_key, table->key_size);
/* memcpy((char *)table->send_entry + table->key_size + 1, (char *)send_data,
*/
/* table->data_size); */
// locking window of target rank with exclusive lock
if (MPI_Win_lock(MPI_LOCK_EXCLUSIVE, dest_rank, 0, table->window) != 0)
return DHT_MPI_ERROR;
for (i = 0; i < table->index_count; i++) {
if (MPI_Get(table->recv_entry, 1 + table->data_size + table->key_size,
MPI_BYTE, dest_rank, table->index[i],
1 + table->data_size + table->key_size, MPI_BYTE,
table->window) != 0)
return DHT_MPI_ERROR;
if (MPI_Win_flush(dest_rank, table->window) != 0)
return DHT_MPI_ERROR;
// increment eviction counter if receiving key doesn't match sending key
// entry has write flag and last index is reached.
if (read_flag(*(char *)table->recv_entry)) {
if (memcmp(send_key, (char *)table->recv_entry + 1, table->key_size) !=
0) {
if (i == (table->index_count) - 1) {
table->evictions += 1;
#ifdef DHT_STATISTICS
table->stats->evictions += 1;
#endif
result = DHT_WRITE_SUCCESS_WITH_EVICTION;
break;
}
} else
break;
} else {
#ifdef DHT_STATISTICS
table->stats->writes_local[dest_rank]++;
#endif
break;
}
}
if (result == DHT_WRITE_SUCCESS_WITH_EVICTION) {
memset((char *)table->send_entry + 1 + table->key_size, '\0',
table->data_size);
} else {
memcpy((char *)table->send_entry + 1 + table->key_size,
(char *)table->recv_entry + 1 + table->key_size, table->data_size);
}
*callback_ret = table->accumulate_callback(
data_size, (char *)send_data, table->data_size,
(char *)table->send_entry + 1 + table->key_size);
// put data to DHT (with last selected index by value i)
if (*callback_ret == 0) {
if (MPI_Put(table->send_entry, 1 + table->data_size + table->key_size,
MPI_BYTE, dest_rank, table->index[i],
1 + table->data_size + table->key_size, MPI_BYTE,
table->window) != 0)
return DHT_MPI_ERROR;
}
// unlock window of target rank
if (MPI_Win_unlock(dest_rank, table->window) != 0)
return DHT_MPI_ERROR;
if (proc) {
*proc = dest_rank;
}
if (index) {
*index = table->index[i];
}
return result;
}
int DHT_write(DHT *table, void *send_key, void *send_data, uint32_t *proc,
uint32_t *index) {
unsigned int dest_rank, i;
int result = DHT_SUCCESS;
@ -146,7 +253,8 @@ int DHT_write(DHT *table, void *send_key, void *send_data) {
1 + table->data_size + table->key_size, MPI_BYTE,
table->window) != 0)
return DHT_MPI_ERROR;
if (MPI_Win_flush(dest_rank, table->window) != 0) return DHT_MPI_ERROR;
if (MPI_Win_flush(dest_rank, table->window) != 0)
return DHT_MPI_ERROR;
// increment eviction counter if receiving key doesn't match sending key
// entry has write flag and last index is reached.
@ -178,12 +286,21 @@ int DHT_write(DHT *table, void *send_key, void *send_data) {
table->window) != 0)
return DHT_MPI_ERROR;
// unlock window of target rank
if (MPI_Win_unlock(dest_rank, table->window) != 0) return DHT_MPI_ERROR;
if (MPI_Win_unlock(dest_rank, table->window) != 0)
return DHT_MPI_ERROR;
if (proc) {
*proc = dest_rank;
}
if (index) {
*index = table->index[i];
}
return result;
}
int DHT_read(DHT *table, void *send_key, void *destination) {
int DHT_read(DHT *table, const void *send_key, void *destination) {
unsigned int dest_rank, i;
#ifdef DHT_STATISTICS
@ -205,7 +322,8 @@ int DHT_read(DHT *table, void *send_key, void *destination) {
1 + table->data_size + table->key_size, MPI_BYTE,
table->window) != 0)
return DHT_MPI_ERROR;
if (MPI_Win_flush(dest_rank, table->window) != 0) return DHT_MPI_ERROR;
if (MPI_Win_flush(dest_rank, table->window) != 0)
return DHT_MPI_ERROR;
// increment read error counter if write flag isn't set ...
if ((read_flag(*(char *)table->recv_entry)) == 0) {
@ -214,7 +332,8 @@ int DHT_read(DHT *table, void *send_key, void *destination) {
table->stats->read_misses += 1;
#endif
// unlock window and return
if (MPI_Win_unlock(dest_rank, table->window) != 0) return DHT_MPI_ERROR;
if (MPI_Win_unlock(dest_rank, table->window) != 0)
return DHT_MPI_ERROR;
return DHT_READ_MISS;
}
@ -227,7 +346,8 @@ int DHT_read(DHT *table, void *send_key, void *destination) {
table->stats->read_misses += 1;
#endif
// unlock window an return
if (MPI_Win_unlock(dest_rank, table->window) != 0) return DHT_MPI_ERROR;
if (MPI_Win_unlock(dest_rank, table->window) != 0)
return DHT_MPI_ERROR;
return DHT_READ_MISS;
}
} else
@ -235,7 +355,8 @@ int DHT_read(DHT *table, void *send_key, void *destination) {
}
// unlock window of target rank
if (MPI_Win_unlock(dest_rank, table->window) != 0) return DHT_MPI_ERROR;
if (MPI_Win_unlock(dest_rank, table->window) != 0)
return DHT_MPI_ERROR;
// if matching key was found copy data into memory of passed pointer
memcpy((char *)destination, (char *)table->recv_entry + table->key_size + 1,
@ -244,6 +365,34 @@ int DHT_read(DHT *table, void *send_key, void *destination) {
return DHT_SUCCESS;
}
int DHT_read_location(DHT *table, uint32_t proc, uint32_t index,
void *destination) {
const uint32_t bucket_size = table->data_size + table->key_size + 1;
#ifdef DHT_STATISTICS
table->stats->r_access++;
#endif
// locking window of target rank with shared lock
if (MPI_Win_lock(MPI_LOCK_SHARED, proc, 0, table->window) != 0)
return DHT_MPI_ERROR;
// receive data
if (MPI_Get(table->recv_entry, bucket_size, MPI_BYTE, proc, index,
bucket_size, MPI_BYTE, table->window) != 0) {
return DHT_MPI_ERROR;
}
// unlock window of target rank
if (MPI_Win_unlock(proc, table->window) != 0)
return DHT_MPI_ERROR;
// if matching key was found copy data into memory of passed pointer
memcpy((char *)destination, (char *)table->recv_entry + 1 + table->key_size,
table->data_size);
return DHT_SUCCESS;
}
int DHT_to_file(DHT *table, const char *filename) {
// open file
MPI_File file;
@ -257,17 +406,15 @@ int DHT_to_file(DHT *table, const char *filename) {
// write header (key_size and data_size)
if (rank == 0) {
if (MPI_File_write(file, &table->key_size, 1, MPI_INT, MPI_STATUS_IGNORE) !=
0)
if (MPI_File_write_shared(file, &table->key_size, 1, MPI_INT,
MPI_STATUS_IGNORE) != 0)
return DHT_FILE_WRITE_ERROR;
if (MPI_File_write(file, &table->data_size, 1, MPI_INT,
MPI_STATUS_IGNORE) != 0)
if (MPI_File_write_shared(file, &table->data_size, 1, MPI_INT,
MPI_STATUS_IGNORE) != 0)
return DHT_FILE_WRITE_ERROR;
}
// seek file pointer behind header for all processes
if (MPI_File_seek_shared(file, DHT_FILEHEADER_SIZE, MPI_SEEK_SET) != 0)
return DHT_FILE_IO_ERROR;
MPI_Barrier(table->communicator);
char *ptr;
int bucket_size = table->key_size + table->data_size + 1;
@ -283,8 +430,12 @@ int DHT_to_file(DHT *table, const char *filename) {
return DHT_FILE_WRITE_ERROR;
}
}
MPI_Barrier(table->communicator);
// close file
if (MPI_File_close(&file) != 0) return DHT_FILE_IO_ERROR;
if (MPI_File_close(&file) != 0)
return DHT_FILE_IO_ERROR;
return DHT_SUCCESS;
}
@ -303,7 +454,8 @@ int DHT_from_file(DHT *table, const char *filename) {
return DHT_FILE_IO_ERROR;
// get file size
if (MPI_File_get_size(file, &f_size) != 0) return DHT_FILE_IO_ERROR;
if (MPI_File_get_size(file, &f_size) != 0)
return DHT_FILE_IO_ERROR;
MPI_Comm_rank(table->communicator, &rank);
@ -322,8 +474,10 @@ int DHT_from_file(DHT *table, const char *filename) {
return DHT_FILE_READ_ERROR;
// compare if written header data and key size matches current sizes
if (*(int *)buffer != table->key_size) return DHT_WRONG_FILE;
if (*(int *)(buffer + 4) != table->data_size) return DHT_WRONG_FILE;
if (*(int *)buffer != table->key_size)
return DHT_WRONG_FILE;
if (*(int *)(buffer + 4) != table->data_size)
return DHT_WRONG_FILE;
// set offset for each process
offset = bucket_size * table->comm_size;
@ -348,14 +502,16 @@ int DHT_from_file(DHT *table, const char *filename) {
// extract key and data and write to DHT
key = buffer;
data = (buffer + table->key_size);
if (DHT_write(table, key, data) == DHT_MPI_ERROR) return DHT_MPI_ERROR;
if (DHT_write(table, key, data, NULL, NULL) == DHT_MPI_ERROR)
return DHT_MPI_ERROR;
// increment current position
cur_pos += offset;
}
free(buffer);
if (MPI_File_close(&file) != 0) return DHT_FILE_IO_ERROR;
if (MPI_File_close(&file) != 0)
return DHT_FILE_IO_ERROR;
return DHT_SUCCESS;
}
@ -377,8 +533,10 @@ int DHT_free(DHT *table, int *eviction_counter, int *readerror_counter) {
return DHT_MPI_ERROR;
*readerror_counter = buf;
}
if (MPI_Win_free(&(table->window)) != 0) return DHT_MPI_ERROR;
if (MPI_Free_mem(table->mem_alloc) != 0) return DHT_MPI_ERROR;
if (MPI_Win_free(&(table->window)) != 0)
return DHT_MPI_ERROR;
if (MPI_Free_mem(table->mem_alloc) != 0)
return DHT_MPI_ERROR;
free(table->recv_entry);
free(table->send_entry);
free(table->index);
@ -407,7 +565,8 @@ int DHT_print_statistics(DHT *table) {
#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
// obtaining all values from all processes in the communicator
if (rank == 0) read_misses = (int *)malloc(table->comm_size * sizeof(int));
if (rank == 0)
read_misses = (int *)malloc(table->comm_size * sizeof(int));
if (MPI_Gather(&table->stats->read_misses, 1, MPI_INT, read_misses, 1,
MPI_INT, 0, table->communicator) != 0)
return DHT_MPI_ERROR;
@ -416,7 +575,8 @@ int DHT_print_statistics(DHT *table) {
return DHT_MPI_ERROR;
table->stats->read_misses = 0;
if (rank == 0) evictions = (int *)malloc(table->comm_size * sizeof(int));
if (rank == 0)
evictions = (int *)malloc(table->comm_size * sizeof(int));
if (MPI_Gather(&table->stats->evictions, 1, MPI_INT, evictions, 1, MPI_INT, 0,
table->communicator) != 0)
return DHT_MPI_ERROR;
@ -425,7 +585,8 @@ int DHT_print_statistics(DHT *table) {
return DHT_MPI_ERROR;
table->stats->evictions = 0;
if (rank == 0) w_access = (int *)malloc(table->comm_size * sizeof(int));
if (rank == 0)
w_access = (int *)malloc(table->comm_size * sizeof(int));
if (MPI_Gather(&table->stats->w_access, 1, MPI_INT, w_access, 1, MPI_INT, 0,
table->communicator) != 0)
return DHT_MPI_ERROR;
@ -434,7 +595,8 @@ int DHT_print_statistics(DHT *table) {
return DHT_MPI_ERROR;
table->stats->w_access = 0;
if (rank == 0) r_access = (int *)malloc(table->comm_size * sizeof(int));
if (rank == 0)
r_access = (int *)malloc(table->comm_size * sizeof(int));
if (MPI_Gather(&table->stats->r_access, 1, MPI_INT, r_access, 1, MPI_INT, 0,
table->communicator) != 0)
return DHT_MPI_ERROR;
@ -443,13 +605,14 @@ int DHT_print_statistics(DHT *table) {
return DHT_MPI_ERROR;
table->stats->r_access = 0;
if (rank == 0) written_buckets = (int *)calloc(table->comm_size, sizeof(int));
if (rank == 0)
written_buckets = (int *)calloc(table->comm_size, sizeof(int));
if (MPI_Reduce(table->stats->writes_local, written_buckets, table->comm_size,
MPI_INT, MPI_SUM, 0, table->communicator) != 0)
return DHT_MPI_ERROR;
if (rank == 0) { // only process with rank 0 will print out results as a
// table
if (rank == 0) { // only process with rank 0 will print out results as a
// table
int sum_written_buckets = 0;
for (int i = 0; i < table->comm_size; i++) {

View File

@ -0,0 +1,300 @@
// Time-stamp: "Last modified 2023-08-01 13:41:57 mluebke"
/*
** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
** Potsdam)
**
** Copyright (C) 2018-2021 Marco De Lucia (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 "poet/DHT_Wrapper.hpp"
#include "poet/DHT_Types.hpp"
#include "poet/HashFunctions.hpp"
#include "poet/Interpolation.hpp"
#include "poet/LookupKey.hpp"
#include "poet/Rounding.hpp"
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <cstring>
#include <iostream>
#include <stdexcept>
#include <vector>
using namespace std;
namespace poet {
DHT_Wrapper::DHT_Wrapper(MPI_Comm dht_comm, std::uint64_t dht_size,
const NamedVector<std::uint32_t> &key_species,
const std::vector<std::int32_t> &key_indices,
uint32_t data_count)
: key_count(key_indices.size()), data_count(data_count),
input_key_elements(key_indices), communicator(dht_comm),
key_species(key_species) {
// initialize DHT object
// key size = count of key elements + timestep
uint32_t key_size = (key_count + 1) * sizeof(Lookup_Keyelement);
uint32_t data_size =
(data_count + input_key_elements.size()) * sizeof(double);
uint32_t buckets_per_process =
static_cast<std::uint32_t>(dht_size / (data_size + key_size));
dht_object = DHT_create(dht_comm, buckets_per_process, data_size, key_size,
&poet::Murmur2_64A);
dht_signif_vector = key_species.getValues();
// this->dht_signif_vector.resize(key_size, DHT_KEY_SIGNIF_DEFAULT);
this->dht_prop_type_vector.resize(key_count, DHT_TYPE_DEFAULT);
const auto key_names = key_species.getNames();
auto tot_h = std::find(key_names.begin(), key_names.end(), "H");
if (tot_h != key_names.end()) {
this->dht_prop_type_vector[tot_h - key_names.begin()] = DHT_TYPE_TOTAL;
}
auto tot_o = std::find(key_names.begin(), key_names.end(), "O");
if (tot_o != key_names.end()) {
this->dht_prop_type_vector[tot_o - key_names.begin()] = DHT_TYPE_TOTAL;
}
auto charge = std::find(key_names.begin(), key_names.end(), "Charge");
if (charge != key_names.end()) {
this->dht_prop_type_vector[charge - key_names.begin()] = DHT_TYPE_CHARGE;
}
}
DHT_Wrapper::~DHT_Wrapper() {
// free DHT
DHT_free(dht_object, NULL, NULL);
}
auto DHT_Wrapper::checkDHT(int length, double dt,
const std::vector<double> &work_package,
std::vector<std::uint32_t> &curr_mapping)
-> const DHT_ResultObject & {
dht_results.length = length;
dht_results.keys.resize(length);
dht_results.results.resize(length);
dht_results.needPhreeqc.resize(length);
std::vector<double> bucket_writer(this->data_count +
input_key_elements.size());
std::vector<std::uint32_t> new_mapping;
// loop over every grid cell contained in work package
for (int i = 0; i < length; i++) {
// point to current grid cell
void *key = (void *)&(work_package[i * this->data_count]);
auto &data = dht_results.results[i];
auto &key_vector = dht_results.keys[i];
// 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
int res =
DHT_read(this->dht_object, key_vector.data(), bucket_writer.data());
switch (res) {
case DHT_SUCCESS:
dht_results.results[i] = inputAndRatesToOutput(bucket_writer);
dht_results.needPhreeqc[i] = false;
this->dht_hits++;
break;
case DHT_READ_MISS:
dht_results.needPhreeqc[i] = true;
new_mapping.push_back(curr_mapping[i]);
dht_results.results[i] = std::vector<double>{
&work_package[i * this->data_count],
&work_package[i * this->data_count] + this->data_count};
// HACK: apply normalization to total H and O in results field of DHT
// dht_results.results[i][0] -= base_totals[0];
// dht_results.results[i][1] -= base_totals[1];
break;
}
}
curr_mapping = std::move(new_mapping);
dht_results.old_values = work_package;
return dht_results;
}
void DHT_Wrapper::fillDHT(int length, const std::vector<double> &work_package) {
// loop over every grid cell contained in work package
dht_results.locations.resize(length);
for (int i = 0; i < length; i++) {
// If true grid cell was simulated, needs to be inserted into dht
if (dht_results.needPhreeqc[i]) {
// check if calcite or dolomite is absent and present, resp.n and vice
// versa in input/output. If this is the case -> Do not write to DHT!
// HACK: hardcoded, should be fixed!
if ((dht_results.old_values[i * this->data_count + 7] == 0) !=
(work_package[i * this->data_count + 7] == 0)) {
dht_results.needPhreeqc[i] = false;
continue;
}
if ((dht_results.old_values[i * this->data_count + 9] == 0) !=
(work_package[i * this->data_count + 9] == 0)) {
dht_results.needPhreeqc[i] = false;
continue;
}
uint32_t proc, index;
const auto &key = dht_results.keys[i];
const auto curr_old_data = std::vector<double>(
dht_results.old_values.begin() + (i * this->data_count),
dht_results.old_values.begin() + ((i + 1) * this->data_count));
const auto curr_new_data = std::vector<double>(
work_package.begin() + (i * this->data_count),
work_package.begin() + ((i + 1) * this->data_count));
const auto data = outputToInputAndRates(curr_old_data, curr_new_data);
// void *data = (void *)&(work_package[i * this->data_count]);
// fuzz data (round, logarithm etc.)
// insert simulated data with fuzzed key into DHT
int res = DHT_write(this->dht_object, (void *)(key.data()),
const_cast<double *>(data.data()), &proc, &index);
dht_results.locations[i] = {proc, index};
// if data was successfully written ...
if ((res != DHT_SUCCESS) && (res == DHT_WRITE_SUCCESS_WITH_EVICTION)) {
dht_evictions++;
}
}
}
}
std::vector<double>
DHT_Wrapper::outputToInputAndRates(const std::vector<double> &old_results,
const std::vector<double> &new_results) {
const int prefix_size = this->input_key_elements.size();
std::vector<double> output(prefix_size + this->data_count);
std::copy(new_results.begin(), new_results.end(),
output.begin() + prefix_size);
for (int i = 0; i < prefix_size; i++) {
const int data_elem_i = input_key_elements[i];
output[i] = old_results[data_elem_i];
output[prefix_size + data_elem_i] -= old_results[data_elem_i];
}
return output;
}
std::vector<double>
DHT_Wrapper::inputAndRatesToOutput(const std::vector<double> &dht_data) {
const int prefix_size = this->input_key_elements.size();
std::vector<double> output{dht_data.begin() + prefix_size, dht_data.end()};
for (int i = 0; i < prefix_size; i++) {
const int data_elem_i = input_key_elements[i];
output[data_elem_i] += dht_data[i];
}
return output;
}
void DHT_Wrapper::resultsToWP(std::vector<double> &work_package) {
for (int i = 0; i < dht_results.length; i++) {
if (!dht_results.needPhreeqc[i]) {
std::copy(dht_results.results[i].begin(), dht_results.results[i].end(),
work_package.begin() + (data_count * i));
}
}
}
int DHT_Wrapper::tableToFile(const char *filename) {
int res = DHT_to_file(dht_object, filename);
return res;
}
int DHT_Wrapper::fileToTable(const char *filename) {
int res = DHT_from_file(dht_object, filename);
if (res != DHT_SUCCESS)
return res;
#ifdef DHT_STATISTICS
DHT_print_statistics(dht_object);
#endif
return DHT_SUCCESS;
}
void DHT_Wrapper::printStatistics() {
int res;
res = DHT_print_statistics(dht_object);
if (res != DHT_SUCCESS) {
// MPI ERROR ... WHAT TO DO NOW?
// RUNNING CIRCLES WHILE SCREAMING
}
}
LookupKey DHT_Wrapper::fuzzForDHT(int var_count, void *key, double dt) {
const auto c_zero_val = std::pow(10, AQUEOUS_EXP);
const Lookup_Keyelement dummy = {.0};
LookupKey vecFuzz(var_count + 1, dummy);
DHT_Rounder rounder;
int totals_i = 0;
// introduce fuzzing to allow more hits in DHT
// loop over every variable of grid cell
for (std::uint32_t i = 0; i < input_key_elements.size(); i++) {
if (input_key_elements[i] == DHT_KEY_INPUT_CUSTOM) {
continue;
}
double curr_key = ((double *)key)[input_key_elements[i]];
if (curr_key != 0) {
if (curr_key < c_zero_val &&
this->dht_prop_type_vector[i] == DHT_TYPE_DEFAULT) {
continue;
}
if (this->dht_prop_type_vector[i] == DHT_TYPE_TOTAL) {
curr_key -= base_totals[totals_i++];
}
vecFuzz[i] =
rounder.round(curr_key, dht_signif_vector[i],
this->dht_prop_type_vector[i] == DHT_TYPE_TOTAL);
}
}
// add timestep to the end of the key as double value
vecFuzz[var_count].fp_element = dt;
return vecFuzz;
}
void poet::DHT_Wrapper::SetSignifVector(std::vector<uint32_t> signif_vec) {
if (signif_vec.size() != this->key_count) {
throw std::runtime_error(
"Significant vector size mismatches count of key elements.");
}
this->dht_signif_vector = signif_vec;
}
} // namespace poet

View File

@ -1,4 +1,4 @@
// Time-stamp: "Last modified 2023-04-24 16:56:18 mluebke"
// Time-stamp: "Last modified 2023-04-24 23:20:55 mluebke"
/*
**-----------------------------------------------------------------------------
** MurmurHash2 was written by Austin Appleby, and is placed in the public

View File

@ -0,0 +1,152 @@
// Time-stamp: "Last modified 2023-08-01 23:18:45 mluebke"
#include "poet/DHT_Wrapper.hpp"
#include "poet/HashFunctions.hpp"
#include "poet/Interpolation.hpp"
#include "poet/LookupKey.hpp"
#include "poet/Rounding.hpp"
#include <algorithm>
#include <array>
#include <cassert>
#include <cstddef>
#include <cstdint>
#include <functional>
#include <iterator>
#include <memory>
#include <mpi.h>
#include <string>
#include <utility>
#include <vector>
extern "C" {
#include "poet/DHT.h"
}
namespace poet {
InterpolationModule::InterpolationModule(
std::uint32_t entries_per_bucket, std::uint64_t size_per_process,
std::uint32_t min_entries_needed, DHT_Wrapper &dht,
const NamedVector<std::uint32_t> &interp_key_signifs,
const std::vector<std::int32_t> &dht_key_indices)
: dht_instance(dht), key_signifs(interp_key_signifs),
key_indices(dht_key_indices), min_entries_needed(min_entries_needed) {
initPHT(this->key_signifs.size(), entries_per_bucket, size_per_process,
dht.getCommunicator());
pht->setSourceDHT(dht.getDHT());
}
void InterpolationModule::initPHT(std::uint32_t key_count,
std::uint32_t entries_per_bucket,
std::uint32_t size_per_process,
MPI_Comm communicator) {
uint32_t key_size = key_count * sizeof(Lookup_Keyelement) + sizeof(double);
uint32_t data_size = sizeof(DHT_Location);
pht = std::make_unique<ProximityHashTable>(
key_size, data_size, entries_per_bucket, size_per_process, communicator);
}
void InterpolationModule::writePairs(const DHT_Wrapper::DHT_ResultObject &in) {
for (int i = 0; i < in.length; i++) {
if (in.needPhreeqc[i]) {
const auto coarse_key = roundKey(in.keys[i]);
pht->writeLocationToPHT(coarse_key, in.locations[i]);
}
}
}
void InterpolationModule::tryInterpolation(
DHT_Wrapper::DHT_ResultObject &dht_results,
std::vector<std::uint32_t> &curr_mapping) {
interp_result.status.resize(dht_results.length, NOT_NEEDED);
interp_result.results.resize(dht_results.length, {});
for (int i = 0; i < dht_results.length; i++) {
if (!dht_results.needPhreeqc[i]) {
interp_result.status[i] = NOT_NEEDED;
continue;
}
auto pht_result =
pht->query(roundKey(dht_results.keys[i]), this->key_signifs.getValues(),
this->min_entries_needed, dht_instance.getInputCount(),
dht_instance.getOutputCount());
int pht_i = 0;
while (pht_i < pht_result.size) {
if (pht_result.size < this->min_entries_needed) {
break;
}
auto in_it = pht_result.in_values.begin() + pht_i;
auto out_it = pht_result.out_values.begin() + pht_i;
bool same_sig_calcite = (pht_result.in_values[pht_i][7] == 0) ==
(dht_results.results[i][7] == 0);
bool same_sig_dolomite = (pht_result.in_values[pht_i][8] == 0) ==
(dht_results.results[i][9] == 0);
if (!same_sig_calcite || !same_sig_dolomite) {
pht_result.size -= 1;
pht_result.in_values.erase(in_it);
pht_result.out_values.erase(out_it);
continue;
}
pht_i += 1;
}
if (pht_result.size < this->min_entries_needed) {
interp_result.status[i] = INSUFFICIENT_DATA;
continue;
}
#ifdef POET_PHT_ADD
this->pht->incrementReadCounter(roundKey(dht_results.keys[i]));
#endif
double start_fc = MPI_Wtime();
// mean water
// double mean_water = 0;
// for (int out_i = 0; out_i < pht_result.size; out_i++) {
// mean_water += pht_result.out_values[out_i][0];
// }
// mean_water /= pht_result.size;
interp_result.results[i] =
f_interpolate(dht_instance.getKeyElements(), dht_results.results[i],
pht_result.in_values, pht_result.out_values);
if (interp_result.results[i][7] < 0 || interp_result.results[i][9] < 0) {
interp_result.status[i] = INSUFFICIENT_DATA;
continue;
}
// interp_result.results[i][0] = mean_water;
this->interp_t += MPI_Wtime() - start_fc;
this->interpolations++;
curr_mapping.erase(std::remove(curr_mapping.begin(), curr_mapping.end(), i),
curr_mapping.end());
dht_results.needPhreeqc[i] = false;
interp_result.status[i] = RES_OK;
}
}
void InterpolationModule::resultsToWP(std::vector<double> &work_package) const {
for (uint32_t i = 0; i < interp_result.status.size(); i++) {
if (interp_result.status[i] == RES_OK) {
const std::size_t length =
interp_result.results[i].end() - interp_result.results[i].begin();
std::copy(interp_result.results[i].begin(),
interp_result.results[i].end(),
work_package.begin() + (length * i));
}
}
}
} // namespace poet

View File

@ -0,0 +1,275 @@
// Time-stamp: "Last modified 2023-08-01 17:11:42 mluebke"
#include "poet/DHT_Wrapper.hpp"
#include "poet/HashFunctions.hpp"
#include "poet/Interpolation.hpp"
#include "poet/LookupKey.hpp"
#include "poet/Rounding.hpp"
#include <cassert>
#include <cstddef>
#include <cstdint>
#include <iostream>
#include <memory>
#include <unordered_set>
#include <vector>
extern "C" {
#include "poet/DHT.h"
}
namespace poet {
ProximityHashTable::ProximityHashTable(uint32_t key_size, uint32_t data_size,
uint32_t entry_count,
uint32_t size_per_process,
MPI_Comm communicator_)
: communicator(communicator_) {
data_size *= entry_count;
data_size += sizeof(bucket_indicator);
#ifdef POET_PHT_ADD
data_size += sizeof(std::uint64_t);
#endif
bucket_store = new char[data_size];
uint32_t buckets_per_process =
static_cast<std::uint32_t>(size_per_process / (data_size + key_size));
this->prox_ht = DHT_create(communicator, buckets_per_process, data_size,
key_size, &poet::Murmur2_64A);
DHT_set_accumulate_callback(this->prox_ht, PHT_callback_function);
}
ProximityHashTable::~ProximityHashTable() {
delete[] bucket_store;
if (prox_ht) {
DHT_free(prox_ht, NULL, NULL);
}
}
int ProximityHashTable::PHT_callback_function(int in_data_size, void *in_data,
int out_data_size,
void *out_data) {
const int max_elements_per_bucket =
static_cast<int>((out_data_size - sizeof(bucket_indicator)
#ifdef POET_PHT_ADD
- sizeof(std::uint64_t)
#endif
) /
in_data_size);
DHT_Location *input = reinterpret_cast<DHT_Location *>(in_data);
bucket_indicator *occupied_buckets =
reinterpret_cast<bucket_indicator *>(out_data);
DHT_Location *pairs = reinterpret_cast<DHT_Location *>(occupied_buckets + 1);
if (*occupied_buckets == max_elements_per_bucket) {
return INTERP_CB_FULL;
}
for (bucket_indicator i = 0; i < *occupied_buckets; i++) {
if (pairs[i] == *input) {
return INTERP_CB_ALREADY_IN;
}
}
pairs[(*occupied_buckets)++] = *input;
return INTERP_CB_OK;
}
void ProximityHashTable::writeLocationToPHT(LookupKey key,
DHT_Location location) {
double start = MPI_Wtime();
// if (localCache[key].first) {
// return;
// }
int ret_val;
int status = DHT_write_accumulate(prox_ht, key.data(), sizeof(location),
&location, NULL, NULL, &ret_val);
if (status == DHT_WRITE_SUCCESS_WITH_EVICTION) {
this->dht_evictions++;
}
// if (ret_val == INTERP_CB_FULL) {
// localCache(key, {});
// }
this->pht_write_t += MPI_Wtime() - start;
}
const ProximityHashTable::PHT_Result &ProximityHashTable::query(
const LookupKey &key, const std::vector<std::uint32_t> &signif,
std::uint32_t min_entries_needed, std::uint32_t input_count,
std::uint32_t output_count) {
double start_r = MPI_Wtime();
const auto cache_ret = localCache[key];
if (cache_ret.first) {
cache_hits++;
return (lookup_results = cache_ret.second);
}
int res = DHT_read(prox_ht, key.data(), bucket_store);
this->pht_read_t += MPI_Wtime() - start_r;
if (res != DHT_SUCCESS) {
this->lookup_results.size = 0;
return lookup_results;
}
auto *bucket_element_count =
reinterpret_cast<bucket_indicator *>(bucket_store);
auto *bucket_elements =
reinterpret_cast<DHT_Location *>(bucket_element_count + 1);
if (*bucket_element_count < min_entries_needed) {
this->lookup_results.size = 0;
return lookup_results;
}
lookup_results.size = *bucket_element_count;
auto locations = std::vector<DHT_Location>(
bucket_elements, bucket_elements + *(bucket_element_count));
lookup_results.in_values.clear();
lookup_results.in_values.reserve(*bucket_element_count);
lookup_results.out_values.clear();
lookup_results.out_values.reserve(*bucket_element_count);
for (const auto &loc : locations) {
double start_g = MPI_Wtime();
DHT_read_location(source_dht, loc.first, loc.second, dht_buffer.data());
this->pht_gather_dht_t += MPI_Wtime() - start_g;
auto *buffer = reinterpret_cast<double *>(dht_buffer.data());
lookup_results.in_values.push_back(
std::vector<double>(buffer, buffer + input_count));
buffer += input_count;
lookup_results.out_values.push_back(
std::vector<double>(buffer, buffer + output_count));
// if (!similarityCheck(check_key, key, signif)) {
// // TODO: original stored location in PHT was overwritten in DHT.
// Need to
// // handle this!
// lookup_results.size--;
// if (lookup_results.size < min_entries_needed) {
// lookup_results.size = 0;
// break;
// }
// continue;
// }
// auto input = convertKeysFromDHT(buffer_start, dht_key_count);
// // remove timestep from the key
// input.pop_back();
// lookup_results.in_keys.push_back(input);
// auto *data = reinterpret_cast<double *>(buffer + dht_key_count);
// lookup_results.out_values.push_back(
// std::vector<double>(data, data + dht_data_count));
}
if (lookup_results.size != 0) {
localCache(key, lookup_results);
}
return lookup_results;
}
inline bool
ProximityHashTable::similarityCheck(const LookupKey &fine,
const LookupKey &coarse,
const std::vector<uint32_t> &signif) {
PHT_Rounder rounder;
for (int i = 0; i < signif.size(); i++) {
if (!(rounder.round(fine[i], signif[i]) == coarse[i])) {
return false;
}
}
return true;
}
inline std::vector<double>
ProximityHashTable::convertKeysFromDHT(Lookup_Keyelement *keys_in,
std::uint32_t key_size) {
std::vector<double> output(key_size);
DHT_Rounder rounder;
for (int i = 0; i < key_size; i++) {
output[i] = rounder.convert(keys_in[i]);
}
return output;
}
void ProximityHashTable::Cache::operator()(const LookupKey &key,
const PHT_Result val) {
const auto elemIt = this->find(key);
if (elemIt == this->end()) {
if (this->free_mem < 0) {
const LookupKey &to_del = this->lru_queue.back();
const auto elem_d = this->find(to_del);
this->free_mem += elem_d->second.getSize();
this->erase(to_del);
this->keyfinder.erase(to_del);
this->lru_queue.pop_back();
}
this->insert({key, val});
this->lru_queue.emplace_front(key);
this->keyfinder[key] = lru_queue.begin();
this->free_mem -= val.getSize();
return;
}
elemIt->second = val;
}
std::pair<bool, ProximityHashTable::PHT_Result>
ProximityHashTable::Cache::operator[](const LookupKey &key) {
const auto elemIt = this->find(key);
if (elemIt == this->end()) {
return {false, {}};
}
this->lru_queue.splice(lru_queue.begin(), lru_queue, this->keyfinder[key]);
return {true, elemIt->second};
}
#ifdef POET_PHT_ADD
static int PHT_increment_counter(int in_data_size, void *in_data,
int out_data_size, void *out_data) {
char *start = reinterpret_cast<char *>(out_data);
std::uint64_t *counter = reinterpret_cast<std::uint64_t *>(
start + out_data_size - sizeof(std::uint64_t));
*counter += 1;
return 0;
}
void ProximityHashTable::incrementReadCounter(const LookupKey &key) {
auto *old_func_ptr = this->prox_ht->accumulate_callback;
DHT_set_accumulate_callback(prox_ht, PHT_increment_counter);
int ret, dummy;
DHT_write_accumulate(prox_ht, key.data(), 0, NULL, NULL, NULL, &ret);
DHT_set_accumulate_callback(prox_ht, old_func_ptr);
}
#endif
} // namespace poet

286
src/ChemistryModule/WorkerFunctions.cpp Normal file → Executable file
View File

@ -1,10 +1,13 @@
// Time-stamp: "Last modified 2023-07-21 17:22:19 mluebke"
// Time-stamp: "Last modified 2023-08-01 17:22:20 mluebke"
#include "IrmResult.h"
#include "poet/ChemistryModule.hpp"
#include "poet/DHT_Wrapper.hpp"
#include "poet/Interpolation.hpp"
#include <algorithm>
#include <bits/stdint-uintn.h>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <iomanip>
@ -15,6 +18,147 @@
#include <string>
#include <vector>
namespace poet {
//std::vector<double>
//inverseDistanceWeighting(const std::vector<std::int32_t> &to_calc,
// const std::vector<double> &from,
// const std::vector<std::vector<double>> &input,
// const std::vector<std::vector<double>> &output) {
// std::vector<double> results = from;
//
// const std::uint32_t buffer_size = input.size() + 1;
// double buffer[buffer_size];
// double from_rescaled;
//
// const std::uint32_t data_set_n = input.size();
// double rescaled[to_calc.size()][data_set_n + 1];
// double weights[data_set_n];
//
// // rescaling over all key elements
// for (int key_comp_i = 0; key_comp_i < to_calc.size(); key_comp_i++) {
// const auto output_comp_i = to_calc[key_comp_i];
//
// // rescale input between 0 and 1
// for (int point_i = 0; point_i < data_set_n; point_i++) {
// rescaled[key_comp_i][point_i] = input[point_i][key_comp_i];
// }
//
// rescaled[key_comp_i][data_set_n] = from[output_comp_i];
//
// const double min = *std::min_element(rescaled[key_comp_i],
// rescaled[key_comp_i] + data_set_n + 1);
// const double max = *std::max_element(rescaled[key_comp_i],
// rescaled[key_comp_i] + data_set_n + 1);
//
// for (int point_i = 0; point_i < data_set_n; point_i++) {
// rescaled[key_comp_i][point_i] =
// ((max - min) != 0
// ? (rescaled[key_comp_i][point_i] - min) / (max - min)
// : 0);
// }
// rescaled[key_comp_i][data_set_n] =
// ((max - min) != 0 ? (from[output_comp_i] - min) / (max - min) : 0);
// }
//
// // calculate distances for each data set
// double inv_sum = 0;
// for (int point_i = 0; point_i < data_set_n; point_i++) {
// double distance = 0;
// for (int key_comp_i = 0; key_comp_i < to_calc.size(); key_comp_i++) {
// distance += std::pow(
// rescaled[key_comp_i][point_i] - rescaled[key_comp_i][data_set_n], 2);
// }
// weights[point_i] = 1 / std::sqrt(distance);
// assert(!std::isnan(weights[point_i]));
// inv_sum += weights[point_i];
// }
//
// assert(!std::isnan(inv_sum));
//
// // actual interpolation
// // bool has_h = false;
// // bool has_o = false;
//
// for (int key_comp_i = 0; key_comp_i < to_calc.size(); key_comp_i++) {
// const auto output_comp_i = to_calc[key_comp_i];
// double key_delta = 0;
//
// // if (interp_i == 0) {
// // has_h = true;
// // }
//
// // if (interp_i == 1) {
// // has_o = true;
// // }
//
// for (int j = 0; j < data_set_n; j++) {
// key_delta += weights[j] * input[j][key_comp_i];
// }
//
// key_delta /= inv_sum;
//
// results[output_comp_i] = from[output_comp_i] + key_delta;
// }
//
// // if (!has_h) {
// // double new_val = 0;
// // for (int j = 0; j < data_set_n; j++) {
// // new_val += weights[j] * output[j][0];
// // }
// // results[0] = new_val / inv_sum;
// // }
//
// // if (!has_h) {
// // double new_val = 0;
// // for (int j = 0; j < data_set_n; j++) {
// // new_val += weights[j] * output[j][1];
// // }
// // results[1] = new_val / inv_sum;
// // }
//
// // for (std::uint32_t i = 0; i < to_calc.size(); i++) {
// // const std::uint32_t interp_i = to_calc[i];
//
// // // rescale input between 0 and 1
// // for (int j = 0; j < input.size(); j++) {
// // buffer[j] = input[j].at(i);
// // }
//
// // buffer[buffer_size - 1] = from[interp_i];
//
// // const double min = *std::min_element(buffer, buffer + buffer_size);
// // const double max = *std::max_element(buffer, buffer + buffer_size);
//
// // for (int j = 0; j < input.size(); j++) {
// // buffer[j] = ((max - min) != 0 ? (buffer[j] - min) / (max - min) : 1);
// // }
// // from_rescaled =
// // ((max - min) != 0 ? (from[interp_i] - min) / (max - min) : 0);
//
// // double inv_sum = 0;
//
// // // calculate distances for each point
// // for (int i = 0; i < input.size(); i++) {
// // const double distance = std::pow(buffer[i] - from_rescaled, 2);
//
// // buffer[i] = distance > 0 ? (1 / std::sqrt(distance)) : 0;
// // inv_sum += buffer[i];
// // }
// // // calculate new values
// // double new_val = 0;
// // for (int i = 0; i < output.size(); i++) {
// // new_val += buffer[i] * output[i][interp_i];
// // }
// // results[interp_i] = new_val / inv_sum;
// // if (std::isnan(results[interp_i])) {
// // std::cout << "nan with new_val = " << output[0][i] << std::endl;
// // }
// // }
//
// return results;
//}
inline std::string get_string(int root, MPI_Comm communicator) {
int count;
MPI_Bcast(&count, 1, MPI_INT, root, communicator);
@ -52,36 +196,6 @@ void poet::ChemistryModule::WorkerLoop() {
initializeField(dummy);
break;
}
case CHEM_DHT_ENABLE: {
bool enable;
ChemBCast(&enable, 1, MPI_CXX_BOOL);
uint32_t size_mb;
ChemBCast(&size_mb, 1, MPI_UINT32_T);
std::vector<std::string> name_dummy;
SetDHTEnabled(enable, size_mb, name_dummy);
break;
}
case CHEM_DHT_SIGNIF_VEC: {
std::vector<uint32_t> input_vec;
SetDHTSignifVector(input_vec);
break;
}
case CHEM_DHT_SNAPS: {
int type;
ChemBCast(&type, 1, MPI_INT);
SetDHTSnaps(type, get_string(0, this->group_comm));
break;
}
case CHEM_DHT_READ_FILE: {
ReadDHTFile(get_string(0, this->group_comm));
break;
}
case CHEM_WORK_LOOP: {
WorkerProcessPkgs(timings, iteration);
break;
@ -96,12 +210,6 @@ void poet::ChemistryModule::WorkerLoop() {
WorkerMetricsToMaster(type);
break;
}
case CHEM_PROGRESSBAR: {
bool enable;
ChemBCast(&enable, 1, MPI_CXX_BOOL);
setProgressBarPrintout(enable);
break;
}
case CHEM_BREAK_MAIN_LOOP: {
WorkerPostSim(iteration);
loop = false;
@ -204,7 +312,9 @@ void poet::ChemistryModule::WorkerDoWork(MPI_Status &probe_status,
dht->checkDHT(local_work_package_size, dt, vecCurrWP, vecMappingWP);
dht_get_end = MPI_Wtime();
// DHT_Results.ResultsToMapping(vecMappingWP);
if (interp_enabled) {
interp->tryInterpolation(dht->getDHTResults(), vecMappingWP);
}
}
phreeqc_time_start = MPI_Wtime();
@ -218,6 +328,9 @@ void poet::ChemistryModule::WorkerDoWork(MPI_Status &probe_status,
if (dht_enabled) {
dht->resultsToWP(vecCurrWP);
if (interp_enabled) {
interp->resultsToWP(vecCurrWP);
}
}
/* send results to master */
@ -231,6 +344,10 @@ void poet::ChemistryModule::WorkerDoWork(MPI_Status &probe_status,
dht->fillDHT(local_work_package_size, vecCurrWP);
dht_fill_end = MPI_Wtime();
if (interp_enabled) {
interp->writePairs(dht->getDHTResults());
}
timings.dht_get += dht_get_end - dht_get_start;
timings.dht_fill += dht_fill_end - dht_fill_start;
}
@ -245,23 +362,44 @@ void poet::ChemistryModule::WorkerPostIter(MPI_Status &prope_status,
MPI_Recv(NULL, 0, MPI_DOUBLE, 0, LOOP_END, this->group_comm,
MPI_STATUS_IGNORE);
if (this->dht_enabled) {
this->dht->printStatistics();
dht_hits.push_back(dht->getHits());
dht_evictions.push_back(dht->getEvictions());
dht->resetCounter();
if (this->dht_snaps_type == DHT_FILES_ITEREND) {
if (this->interp_enabled) {
interp_calls.push_back(interp->getInterpolationCount());
interp->resetCounter();
interp->writePHTStats();
}
if (this->dht_snaps_type == DHT_SNAPS_ITEREND) {
WorkerWriteDHTDump(iteration);
if (this->interp_enabled) {
std::stringstream out;
out << this->dht_file_out_dir << "/iter_" << std::setfill('0')
<< std::setw(this->file_pad) << iteration << ".pht";
interp->dumpPHTState(out.str());
}
}
}
}
void poet::ChemistryModule::WorkerPostSim(uint32_t iteration) {
if (this->dht_enabled && this->dht_snaps_type == DHT_FILES_SIMEND) {
if (this->dht_enabled && this->dht_snaps_type == DHT_SNAPS_SIMEND) {
WorkerWriteDHTDump(iteration);
if (this->interp_enabled) {
std::stringstream out;
out << this->dht_file_out_dir << "/iter_" << std::setfill('0')
<< std::setw(this->file_pad) << iteration << ".pht";
interp->dumpPHTState(out.str());
}
}
}
void poet::ChemistryModule::WorkerWriteDHTDump(uint32_t iteration) {
std::stringstream out;
out << this->dht_file_out_dir << "/iter_" << std::setfill('0') << std::setw(3)
<< iteration << ".dht";
out << this->dht_file_out_dir << "/iter_" << std::setfill('0')
<< std::setw(this->file_pad) << iteration << ".dht";
int res = dht->tableToFile(out.str().c_str());
if (res != DHT_SUCCESS && this->comm_rank == 2)
std::cerr
@ -348,6 +486,26 @@ void poet::ChemistryModule::WorkerPerfToMaster(int type,
this->group_comm);
break;
}
case WORKER_IP_WRITE: {
double val = interp->getPHTWriteTime();
MPI_Gather(&val, 1, MPI_DOUBLE, NULL, 1, MPI_DOUBLE, 0, this->group_comm);
break;
}
case WORKER_IP_READ: {
double val = interp->getPHTReadTime();
MPI_Gather(&val, 1, MPI_DOUBLE, NULL, 1, MPI_DOUBLE, 0, this->group_comm);
break;
}
case WORKER_IP_GATHER: {
double val = interp->getDHTGatherTime();
MPI_Gather(&val, 1, MPI_DOUBLE, NULL, 1, MPI_DOUBLE, 0, this->group_comm);
break;
}
case WORKER_IP_FC: {
double val = interp->getInterpolationTime();
MPI_Gather(&val, 1, MPI_DOUBLE, NULL, 1, MPI_DOUBLE, 0, this->group_comm);
break;
}
default: {
throw std::runtime_error("Unknown perf type in master's message.");
}
@ -355,24 +513,46 @@ void poet::ChemistryModule::WorkerPerfToMaster(int type,
}
void poet::ChemistryModule::WorkerMetricsToMaster(int type) {
uint32_t value;
MPI_Comm worker_comm = dht->getCommunicator();
int worker_rank;
MPI_Comm_rank(worker_comm, &worker_rank);
MPI_Comm &group_comm = this->group_comm;
auto reduce_and_send = [&worker_rank, &worker_comm, &group_comm](
std::vector<std::uint32_t> &send_buffer, int tag) {
std::vector<uint32_t> to_master(send_buffer.size());
MPI_Reduce(send_buffer.data(), to_master.data(), send_buffer.size(),
MPI_UINT32_T, MPI_SUM, 0, worker_comm);
if (worker_rank == 0) {
MPI_Send(to_master.data(), to_master.size(), MPI_UINT32_T, 0, tag,
group_comm);
}
};
switch (type) {
case WORKER_DHT_HITS: {
value = dht->getHits();
break;
}
case WORKER_DHT_MISS: {
value = dht->getMisses();
reduce_and_send(dht_hits, WORKER_DHT_HITS);
break;
}
case WORKER_DHT_EVICTIONS: {
value = dht->getEvictions();
reduce_and_send(dht_evictions, WORKER_DHT_EVICTIONS);
break;
}
case WORKER_IP_CALLS: {
reduce_and_send(interp_calls, WORKER_IP_CALLS);
return;
}
case WORKER_PHT_CACHE_HITS: {
std::vector<std::uint32_t> input = this->interp->getPHTLocalCacheHits();
reduce_and_send(input, WORKER_PHT_CACHE_HITS);
return;
}
default: {
throw std::runtime_error("Unknown perf type in master's message.");
}
}
MPI_Gather(&value, 1, MPI_UINT32_T, NULL, 1, MPI_UINT32_T, 0,
this->group_comm);
}
} // namespace poet

0
src/DataStructures/Field.cpp Normal file → Executable file
View File

0
src/DiffusionModule.cpp Normal file → Executable file
View File

0
src/Grid.cpp Normal file → Executable file
View File

73
src/SimParams.cpp Normal file → Executable file
View File

@ -56,14 +56,6 @@ poet::GridParams::s_GridParams(RInside &R) {
(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"));
this->init_df =
Rcpp::as<Rcpp::DataFrame>(R.parseEval("mysetup$grid$init_cell"));
this->props =
Rcpp::as<std::vector<std::string>>(R.parseEval("mysetup$grid$props"));
this->input_script =
Rcpp::as<std::string>(R.parseEval("mysetup$grid$input_script"));
this->database_path =
Rcpp::as<std::string>(R.parseEval("mysetup$grid$database"));
}
poet::DiffusionParams::s_DiffusionParams(RInside &R) {
@ -82,20 +74,20 @@ poet::DiffusionParams::s_DiffusionParams(RInside &R) {
Rcpp::as<Rcpp::DataFrame>(R.parseEval("mysetup$diffusion$vecinj_index"));
}
poet::ChemistryParams::s_ChemistryParams(RInside &R) {
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 (Rcpp::as<bool>(
R.parseEval("'dht_species' %in% names(mysetup$chemistry)"))) {
this->dht_species = Rcpp::as<std::vector<std::string>>(
R.parseEval("mysetup$chemistry$dht_species"));
if (R.checkIfExists("dht_species", "mysetup$chemistry")) {
this->dht_signifs =
R.wrapNamedVector<std::uint32_t>("mysetup$chemistry$dht_species");
}
if (Rcpp::as<bool>(
R.parseEval("'dht_signif' %in% names(mysetup$chemistry)"))) {
this->dht_signif = Rcpp::as<std::vector<std::uint32_t>>(
R.parseEval("mysetup$chemistry$dht_signif"));
if (R.checkIfExists("pht_species", "mysetup$chemistry")) {
this->pht_signifs =
R.wrapNamedVector<std::uint32_t>("mysetup$chemistry$pht_species");
}
}
@ -104,7 +96,7 @@ SimParams::SimParams(int world_rank_, int world_size_) {
this->simparams.world_size = world_size_;
}
int SimParams::parseFromCmdl(char *argv[], RInside &R) {
int SimParams::parseFromCmdl(char *argv[], RInsidePOET &R) {
// initialize argh object
argh::parser cmdl(argv);
@ -144,32 +136,26 @@ int SimParams::parseFromCmdl(char *argv[], RInside &R) {
simparams.print_progressbar = cmdl[{"P", "progress"}];
/*Parse DHT arguments*/
simparams.dht_enabled = cmdl["dht"];
chem_params.use_dht = cmdl["dht"];
// cout << "CPP: DHT is " << ( dht_enabled ? "ON" : "OFF" ) << '\n';
if (simparams.dht_enabled) {
cmdl("dht-strategy", 0) >> simparams.dht_strategy;
// cout << "CPP: DHT strategy is " << dht_strategy << endl;
cmdl("dht-signif", 5) >> simparams.dht_significant_digits;
// cout << "CPP: DHT significant digits = " << dht_significant_digits <<
// endl;
simparams.dht_log = !(cmdl["dht-nolog"]);
// cout << "CPP: DHT logarithm before rounding: " << ( dht_logarithm ? "ON"
// : "OFF" ) << endl;
cmdl("dht-size", DHT_SIZE_PER_PROCESS_MB) >> simparams.dht_size_per_process;
if (chem_params.use_dht) {
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) >> simparams.dht_snaps;
cmdl("dht-snaps", 0) >> chem_params.dht_snaps;
cmdl("dht-file") >> dht_file;
cmdl("dht-file") >> chem_params.dht_file;
}
/*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"];
@ -177,25 +163,27 @@ int SimParams::parseFromCmdl(char *argv[], RInside &R) {
cout << "CPP: Complete results storage is "
<< (simparams.store_result ? "ON" : "OFF") << endl;
cout << "CPP: Work Package Size: " << simparams.wp_size << endl;
cout << "CPP: DHT is " << (simparams.dht_enabled ? "ON" : "OFF") << '\n';
cout << "CPP: DHT is " << (chem_params.use_dht ? "ON" : "OFF") << '\n';
if (simparams.dht_enabled) {
if (chem_params.use_dht) {
cout << "CPP: DHT strategy is " << simparams.dht_strategy << endl;
cout << "CPP: DHT key default digits (ignored if 'signif_vector' is "
"defined) = "
<< simparams.dht_significant_digits << endl;
cout << "CPP: DHT logarithm before rounding: "
<< (simparams.dht_log ? "ON" : "OFF") << endl;
cout << "CPP: DHT size per process (Byte) = "
<< simparams.dht_size_per_process << endl;
cout << "CPP: DHT save snapshots is " << simparams.dht_snaps << endl;
cout << "CPP: DHT load file is " << dht_file << endl;
cout << "CPP: DHT size per process (Byte) = " << chem_params.dht_size
<< endl;
cout << "CPP: DHT save snapshots is " << chem_params.dht_snaps << endl;
cout << "CPP: DHT load file is " << chem_params.dht_file << endl;
}
}
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;
@ -210,12 +198,15 @@ int SimParams::parseFromCmdl(char *argv[], RInside &R) {
// work package size
R["work_package_size"] = simparams.wp_size;
// dht enabled?
R["dht_enabled"] = simparams.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;
}

0
test/CMakeLists.txt Normal file → Executable file
View File

Some files were not shown because too many files have changed in this diff Show More