diff --git a/.devcontainer/Dockerfile b/.devcontainer/Dockerfile index 8016eb576..f3ad8d911 100644 --- a/.devcontainer/Dockerfile +++ b/.devcontainer/Dockerfile @@ -2,11 +2,19 @@ FROM mcr.microsoft.com/vscode/devcontainers/base:debian RUN sudo apt-get update && export DEBIAN_FRONTEND=noninteractive \ && sudo apt-get install -y \ - cmake-curses-gui \ - clangd \ + cmake \ git \ libeigen3-dev \ libopenmpi-dev \ - r-cran-rcpp \ - r-cran-rinside + r-base-dev +RUN git clone https://github.com/doctest/doctest.git /doctest \ + && cd /doctest \ + && mkdir build \ + && cd build \ + && cmake .. \ + && make install \ + && cd / \ + && rm -rf /doctest + +RUN /usr/bin/R -q -e "install.packages(c('Rcpp', 'RInside'))" \ No newline at end of file diff --git a/.gitignore b/.gitignore index 1e597b203..0e1d0720f 100644 --- a/.gitignore +++ b/.gitignore @@ -142,4 +142,5 @@ vignettes/*.pdf build/ /.cache/ -.vscode \ No newline at end of file +.vscode +.codechecker \ No newline at end of file diff --git a/.gitmodules b/.gitmodules index 744dfdb41..26a3acc7f 100644 --- a/.gitmodules +++ b/.gitmodules @@ -2,9 +2,6 @@ path = ext/tug url = ../tug.git -[submodule "ext/phreeqcrm"] - path = ext/phreeqcrm - url = ../phreeqcrm-gfz.git -[submodule "ext/doctest"] - path = ext/doctest - url = https://github.com/doctest/doctest.git +[submodule "ext/iphreeqc"] + path = ext/iphreeqc + url = ../iphreeqc.git diff --git a/CMakeLists.txt b/CMakeLists.txt index 8faf514e6..f2f4434bf 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,12 +1,11 @@ -# Version 3.9+ offers new MPI package variables -cmake_minimum_required(VERSION 3.9) +cmake_minimum_required(VERSION 3.14) project(POET LANGUAGES CXX C DESCRIPTION "A coupled reactive transport simulator") # specify the C++ standard -set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD 20) set(CMAKE_CXX_STANDARD_REQUIRED True) set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE) @@ -16,8 +15,6 @@ list(APPEND CMAKE_MODULE_PATH "${POET_SOURCE_DIR}/CMake") get_poet_version() -# set(GCC_CXX_FLAGS "-D STRICT_R_HEADERS") add_definitions(${GCC_CXX_FLAGS}) - find_package(MPI REQUIRED) find_package(RRuntime REQUIRED) @@ -29,15 +26,14 @@ add_subdirectory(bench) set(TUG_ENABLE_TESTING OFF CACHE BOOL "" FORCE) add_subdirectory(ext/tug EXCLUDE_FROM_ALL) -add_subdirectory(ext/phreeqcrm EXCLUDE_FROM_ALL) +add_subdirectory(ext/iphreeqc EXCLUDE_FROM_ALL) option(POET_ENABLE_TESTING "Build test suite for POET" OFF) if (POET_ENABLE_TESTING) - add_subdirectory(ext/doctest EXCLUDE_FROM_ALL) add_subdirectory(test) endif() option(BUILD_DOC "Build documentation with doxygen" OFF) -add_subdirectory(docs) +add_subdirectory(docs) \ No newline at end of file diff --git a/README.md b/README.md index a273f7f10..e6535ed5f 100644 --- a/README.md +++ b/README.md @@ -20,10 +20,10 @@ pages](https://naaice.git-pages.gfz-potsdam.de/poet). The following external header library is shipped with POET: - **argh** - https://github.com/adishavit/argh (BSD license) -- **PhreeqcRM** with patches from GFZ - - https://www.usgs.gov/software/phreeqc-version-3 - - https://git.gfz-potsdam.de/mluebke/phreeqcrm-gfz -- **tug** - https://git.gfz-potsdam.de/sec34/tug +- **IPhreeqc** with patches from GFZ - + https://github.com/usgs-coupled/iphreeqc - + https://git.gfz-potsdam.de/naaice/iphreeqc +- **tug** - https://git.gfz-potsdam.de/naaice/tug ## Installation @@ -35,6 +35,7 @@ To compile POET you need several software to be installed: - MPI-Implementation (tested with OpenMPI and MVAPICH) - R language and environment - CMake 3.9+ +- Eigen3 3.4+ (required by `tug`) - *optional*: `doxygen` with `dot` bindings for documentiation The following R libraries must then be installed, which will get the @@ -107,58 +108,50 @@ The correspondending directory tree would look like this: ```sh poet ├── bin -│ └── poet -├── R_lib -│ └── kin_r_library.R +│   ├── poet +│   └── poet_init └── share └── poet - └── bench - ├── barite - │ ├── barite_interp_eval.R - │ ├── barite.pqi - │ ├── barite.R - │ └── db_barite.dat - ├── dolo - │ ├── dolo_diffu_inner_large.R - │ ├── dolo_diffu_inner.R - │ ├── dolo_inner.pqi - │ ├── dolo_interp_long.R - │ └── phreeqc_kin.dat - └── surfex - ├── ExBase.pqi - ├── ex.R - ├── SMILE_2021_11_01_TH.dat - ├── SurfExBase.pqi - └── surfex.R + ├── barite + │   ├── barite_200.rds + │   ├── barite_200_rt.R + │   ├── barite_het.rds + │   └── barite_het_rt.R + ├── dolo + │   ├── dolo_inner_large.rds + │   ├── dolo_inner_large_rt.R + │   ├── dolo_interp.rds + │   └── dolo_interp_rt.R + └── surfex + ├── PoetEGU_surfex_500.rds + └── PoetEGU_surfex_500_rt.R ``` -The R libraries will be loaded at runtime and the paths are hardcoded -absolute paths inside `poet.cpp`. So, if you consider to move -`bin/poet` either change paths of the R source files and recompile -POET or also move `R_lib/*` relative to the binary. +With the installation of POET, two executables are provided: + - `poet` - the main executable to run simulations + - `poet_init` - a preprocessor to generate input files for POET from R scripts -The benchmarks consist of input scripts, which are provided as .R files. -Additionally, Phreeqc scripts and their corresponding databases are required, -stored as .pqi and .dat files, respectively. +Preprocessed benchmarks can be found in the `share/poet` directory with an +according *runtime* setup. More on those files and how to create them later. ## Running -Run POET by `mpirun ./poet ` +Run POET by `mpirun ./poet [OPTIONS] ` where: -- **OPTIONS** - runtime parameters (explained below) -- **SIMFILE** - simulation described as R script (e.g. - `/share/poet/bench/dolo/dolo_interp_long.R`) +- **OPTIONS** - POET options (explained below) +- **RUNFILE** - Runtime parameters described as R script +- **SIMFILE** - Simulation input prepared by `poet_init` - **OUTPUT_DIRECTORY** - path, where all output of POET should be stored -### Runtime options +### POET options The following parameters can be set: | Option | Value | Description | |-----------------------------|--------------|--------------------------------------------------------------------------------------------------------------------------| | **--work-package-size=** | _1..n_ | size of work packages (defaults to _5_) | -| **--ignore-result** | | disables store of simulation resuls | +| **-P, --progress** | | show progress bar | | **--dht** | | enabling DHT usage (defaults to _OFF_) | | **--dht-strategy=** | _0-1_ | change DHT strategy. **NOT IMPLEMENTED YET** (Defaults to _0_) | | **--dht-size=** | _1-n_ | size of DHT per process involved in megabyte (defaults to _1000 MByte_) | @@ -180,14 +173,16 @@ Following values can be set: ### Example: Running from scratch -We will continue the above example and start a simulation with -`dolo_diffu_inner.R`. As transport a simple fixed-coefficient diffusion is used. -It's a 2D, 100x100 grid, simulating 10 time steps. To start the simulation with -4 processes `cd` into your previously installed POET-dir -`/bin` and run: +We will continue the above example and start a simulation with *barite_het*, +which simulation files can be found in +`/share/poet/barite/barite_het*`. As transport a heterogeneous +diffusion is used. It's a small 2D grid, 2x5 grid, simulating 50 time steps with +a time step size of 100 seconds. To start the simulation with 4 processes `cd` +into your previously installed POET-dir `/bin` and run: ```sh -mpirun -n 4 ./poet ../share/poet/bench/dolo/dolo_diffu_inner.R/ output +cp ../share/poet/barite/barite_het* . +mpirun -n 4 ./poet barite_het_rt.R barite_het.rds output ``` After a finished simulation all data generated by POET will be found @@ -200,9 +195,32 @@ produced. This is done by appending the `--dht-snaps=` option. The resulting call would look like this: ```sh -mpirun -n 4 ./poet --dht --dht-snaps=2 ../share/poet/bench/dolo/dolo_diffu_inner.R/ output +mpirun -n 4 ./poet --dht --dht-snaps=2 barite_het_rt.R barite_het.rds output ``` +## Defining a model + +In order to provide a model to POET, you need to setup a R script which can then +be used by `poet_init` to generate the simulation input. Which parameters are +required can be found in the +[Wiki](https://git.gfz-potsdam.de/naaice/poet/-/wikis/Initialization). We try to +keep the document up-to-date. However, if you encounter missing information or +need help, please get in touch with us via the issue tracker or E-Mail. + +`poet_init` can be used as follows: + +```sh +./poet_init [-o, --output output_file] [-s, --setwd] +``` + +where: + +- **output** - name of the output file (defaults to the input file name + with the extension `.rds`) +- **setwd** - set the working directory to the directory of the input file (e.g. + to allow relative paths in the input script). However, the output file + will be stored in the directory from which `poet_init` was called. + ## About the usage of MPI_Wtime() Implemented time measurement functions uses `MPI_Wtime()`. Some diff --git a/R_lib/init_r_lib.R b/R_lib/init_r_lib.R new file mode 100644 index 000000000..640c0e07b --- /dev/null +++ b/R_lib/init_r_lib.R @@ -0,0 +1,56 @@ +pqc_to_grid <- function(pqc_in, grid) { + # Convert the input DataFrame to a matrix + pqc_in <- as.matrix(pqc_in) + + # Flatten the matrix into a vector + id_vector <- as.numeric(t(grid)) + + # Find the matching rows in the matrix + row_indices <- match(id_vector, pqc_in[, "ID"]) + + # Extract the matching rows from pqc_in to size of grid matrix + result_mat <- pqc_in[row_indices, ] + + # Convert the result matrix to a data frame + res_df <- as.data.frame(result_mat) + + # Remove all columns which only contain NaN + res_df <- res_df[, colSums(is.na(res_df)) != nrow(res_df)] + + # Remove row names + rownames(res_df) <- NULL + + return(res_df) +} + +resolve_pqc_bound <- function(pqc_mat, transport_spec, id) { + df <- as.data.frame(pqc_mat, check.names = FALSE) + value <- df[df$ID == id, transport_spec] + + if (is.nan(value)) { + value <- 0 + } + + return(value) +} + +add_missing_transport_species <- function(init_grid, new_names) { + # add 'ID' to new_names front, as it is not a transport species but required + new_names <- c("ID", new_names) + sol_length <- length(new_names) + + new_grid <- data.frame(matrix(0, nrow = nrow(init_grid), ncol = sol_length)) + names(new_grid) <- new_names + + matching_cols <- intersect(names(init_grid), new_names) + + # Copy matching columns from init_grid to new_grid + new_grid[, matching_cols] <- init_grid[, matching_cols] + + + # Add missing columns to new_grid + append_df <- init_grid[, !(names(init_grid) %in% new_names)] + new_grid <- cbind(new_grid, append_df) + + return(new_grid) +} \ No newline at end of file diff --git a/R_lib/kin_r_library.R b/R_lib/kin_r_library.R index 1e7f803ed..58941a032 100644 --- a/R_lib/kin_r_library.R +++ b/R_lib/kin_r_library.R @@ -15,239 +15,105 @@ ### this program; if not, write to the Free Software Foundation, Inc., 51 ### Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. - -## Simple function to check file extension. It is needed to check if -## the GridFile is SUM (MUFITS format) or rds/RData -FileExt <- function(x) { - pos <- regexpr("\\.([[:alnum:]]+)$", x) - ifelse(pos > -1L, substring(x, pos + 1L), "") -} - -master_init <- function(setup) { - msgm("Process with rank 0 reading GRID properties") - +master_init <- function(setup, out_dir, init_field) { ## Setup the directory where we will store the results - verb <- FALSE - if (local_rank == 0) { - verb <- TRUE ## verbosity loading MUFITS results - if (!dir.exists(fileout)) { - dir.create(fileout) - msgm("created directory ", fileout) - } else { - msgm("dir ", fileout, " already exists, I will overwrite!") - } - if (!exists("store_result")) { - msgm("store_result doesn't exist!") - } else { - msgm("store_result is ", store_result) - } + if (!dir.exists(out_dir)) { + dir.create(out_dir) + msgm("created directory ", out_dir) } else { - + msgm("dir ", out_dir, " already exists, I will overwrite!") + } + if (is.null(setup$store_result)) { + msgm("store_result doesn't exist!") + } else { + msgm("store_result is ", setup$store_result) } setup$iter <- 1 - setup$maxiter <- setup$iterations setup$timesteps <- setup$timesteps + setup$maxiter <- length(setup$timesteps) + setup$iterations <- setup$maxiter setup$simulation_time <- 0 if (is.null(setup[["store_result"]])) { setup$store_result <- TRUE } - + if (setup$store_result) { + init_field_out <- paste0(out_dir, "/iter_0.rds") + init_field <- data.frame(init_field, check.names = FALSE) + saveRDS(init_field, file = init_field_out) + msgm("Stored initial field in ", init_field_out) if (is.null(setup[["out_save"]])) { setup$out_save <- seq(1, setup$iterations) } } - + + setup$out_dir <- out_dir + return(setup) } ## This function, called only by master, stores on disk the last ## calculated time step if store_result is TRUE and increments the ## iteration counter -master_iteration_end <- function(setup) { +master_iteration_end <- function(setup, state_T, state_C) { iter <- setup$iter + # print(iter) ## max digits for iterations - dgts <- as.integer(ceiling(log10(setup$iterations + 1))) - ## string format to use in sprintf + dgts <- as.integer(ceiling(log10(setup$maxiter))) + ## string format to use in sprintf fmt <- paste0("%0", dgts, "d") - + ## Write on disk state_T and state_C after every iteration ## comprised in setup$out_save if (setup$store_result) { if (iter %in% setup$out_save) { - nameout <- paste0(fileout, "/iter_", sprintf(fmt=fmt, iter), ".rds") - info <- list( - tr_req_dt = as.integer(setup$req_dt) - ## tr_allow_dt = setup$allowed_dt, - ## tr_inniter = as.integer(setup$inniter) - ) + nameout <- paste0(setup$out_dir, "/iter_", sprintf(fmt = fmt, iter), ".rds") + + state_T <- data.frame(state_T, check.names = FALSE) + state_C <- data.frame(state_C, check.names = FALSE) + saveRDS(list( - T = setup$state_T, C = setup$state_C, - simtime = as.integer(setup$simtime), - tr_info = info + T = state_T, + C = state_C, + simtime = as.integer(setup$simulation_time) ), file = nameout) msgm("results stored in <", nameout, ">") } } - msgm("done iteration", iter, "/", setup$maxiter) + ## Add last time step to simulation time + setup$simulation_time <- setup$simulation_time + setup$timesteps[iter] + + msgm("done iteration", iter, "/", length(setup$timesteps)) setup$iter <- setup$iter + 1 return(setup) } -## function for the workers to compute chemistry through PHREEQC -slave_chemistry <- function(setup, data) { - base <- setup$base - first <- setup$first - prop <- setup$prop - immobile <- setup$immobile - kin <- setup$kin - ann <- setup$ann - - iter <- setup$iter - timesteps <- setup$timesteps - dt <- timesteps[iter] - - state_T <- data ## not the global field, but the work-package - - ## treat special H+/pH, e-/pe cases - state_T <- RedModRphree::Act2pH(state_T) - - ## reduction of the problem - if (setup$reduce) { - reduced <- ReduceStateOmit(state_T, omit = setup$ann) - } else { - reduced <- state_T - } - - ## form the PHREEQC input script for the current work package - inplist <- SplitMultiKin( - data = reduced, procs = 1, base = base, first = first, - ann = ann, prop = prop, minerals = immobile, kin = kin, dt = dt - ) - - ## if (local_rank==1 & iter==1) - ## RPhreeWriteInp("FirstInp", inplist) - - tmpC <- RunPQC(inplist, procs = 1, second = TRUE) - - ## recompose after the reduction - if (setup$reduce) { - state_C <- RecomposeState(tmpC, reduced) - } else { - state_C <- tmpC - } - - ## the next line is needed since we don't need all columns of - ## PHREEQC output - return(state_C[, prop]) -} - -## This function, called by master -master_chemistry <- function(setup, data) { - state_T <- setup$state_T - - msgm(" chemistry iteration", setup$iter) - - ## treat special H+/pH, e-/pe cases - state_T <- RedModRphree::Act2pH(state_T) - - ## reduction of the problem - if (setup$reduce) { - reduced <- ReduceStateOmit(state_T, omit = setup$ann) - } else { - reduced <- state_T - } - - ## inject data from workers - res_C <- data - - rownames(res_C) <- NULL - - ## print(res_C) - - if (nrow(res_C) > nrow(reduced)) { - res_C <- res_C[seq(2, nrow(res_C), by = 2), ] - } - - ## recompose after the reduction - if (setup$reduce) { - state_C <- RecomposeState(res_C, reduced) - } else { - state_C <- res_C - } - - setup$state_C <- state_C - setup$reduced <- reduced - - return(setup) -} - - -## Adapted version for "reduction" -ReduceStateOmit <- function(data, omit = NULL, sign = 6) { - require(mgcv) - - rem <- colnames(data) - if (is.list(omit)) { - indomi <- match(names(omit), colnames(data)) - datao <- data[, -indomi] - } else { - datao <- data - } - - datao <- signif(datao, sign) - red <- mgcv::uniquecombs(datao) - inds <- attr(red, "index") - now <- ncol(red) - - - ## reattach the omitted column(s) - ## FIXME: control if more than one ann is present - if (is.list(omit)) { - red <- cbind(red, rep(data[1, indomi], nrow(red))) - - colnames(red)[now + 1] <- names(omit) - - ret <- red[, colnames(data)] - } else { - ret <- red - } - rownames(ret) <- NULL - attr(ret, "index") <- inds - return(ret) -} - - ## Attach the name of the calling function to the message displayed on ## R's stdout msgm <- function(...) { - if (local_rank == 0) { - fname <- as.list(sys.call(-1))[[1]] - prefix <- paste0("R: ", fname, " ::") - cat(paste(prefix, ..., "\n")) - } + prefix <- paste0("R: ") + cat(paste(prefix, ..., "\n")) invisible() } ## Function called by master R process to store on disk all relevant ## parameters for the simulation -StoreSetup <- function(setup) { - +StoreSetup <- function(setup, filesim, out_dir) { to_store <- vector(mode = "list", length = 4) ## names(to_store) <- c("Sim", "Flow", "Transport", "Chemistry", "DHT") names(to_store) <- c("Sim", "Transport", "DHT", "Cmdline") - + ## read the setup R file, which is sourced in kin.cpp tmpbuff <- file(filesim, "r") setupfile <- readLines(tmpbuff) close.connection(tmpbuff) - + to_store$Sim <- setupfile - + ## to_store$Flow <- list( ## snapshots = setup$snapshots, ## gridfile = setup$gridfile, @@ -282,22 +148,22 @@ StoreSetup <- function(setup) { to_store$DHT <- FALSE } - if (dht_enabled) { - to_store$DHT <- list( - enabled = dht_enabled, - log = dht_log - #signif = dht_final_signif, - #proptype = dht_final_proptype - ) - } else { - to_store$DHT <- FALSE - } + if (dht_enabled) { + to_store$DHT <- list( + enabled = dht_enabled, + log = dht_log + # signif = dht_final_signif, + # proptype = dht_final_proptype + ) + } else { + to_store$DHT <- FALSE + } - saveRDS(to_store, file = paste0(fileout, "/setup.rds")) - msgm("initialization stored in ", paste0(fileout, "/setup.rds")) + saveRDS(to_store, file = paste0(fileout, "/setup.rds")) + msgm("initialization stored in ", paste0(fileout, "/setup.rds")) } GetWorkPackageSizesVector <- function(n_packages, package_size, len) { - ids <- rep(1:n_packages, times=package_size, each = 1)[1:len] + ids <- rep(1:n_packages, times = package_size, each = 1)[1:len] return(as.integer(table(ids))) } diff --git a/apps/CMakeLists.txt b/apps/CMakeLists.txt new file mode 100644 index 000000000..928a20444 --- /dev/null +++ b/apps/CMakeLists.txt @@ -0,0 +1,4 @@ +file(GLOB INIT_SRCS CONFIGURE_DEPENDS "initializer/*.cpp") + +add_executable(poet_initializer ${INIT_SRCS}) +target_link_libraries(poet_initializer RRuntime tug) \ No newline at end of file diff --git a/apps/initializer/main.cpp b/apps/initializer/main.cpp new file mode 100644 index 000000000..5297158c5 --- /dev/null +++ b/apps/initializer/main.cpp @@ -0,0 +1,3 @@ +#include + +int main(int argc, char **argv) {} \ No newline at end of file diff --git a/bench/CMakeLists.txt b/bench/CMakeLists.txt index 6b5e9b9f3..17045c31b 100644 --- a/bench/CMakeLists.txt +++ b/bench/CMakeLists.txt @@ -1,3 +1,44 @@ + +function(ADD_BENCH_TARGET TARGET POET_BENCH_LIST RT_FILES OUT_PATH) + set(bench_install_dir share/poet/${OUT_PATH}) + + # create empty list + set(OUT_FILES_LIST "") + + foreach(BENCH_FILE ${${POET_BENCH_LIST}}) + get_filename_component(BENCH_NAME ${BENCH_FILE} NAME_WE) + set(OUT_FILE ${CMAKE_CURRENT_BINARY_DIR}/${BENCH_NAME}.rds) + + add_custom_command( + OUTPUT ${OUT_FILE} + COMMAND $ -o ${OUT_FILE} -s ${CMAKE_CURRENT_SOURCE_DIR}/${BENCH_FILE} + COMMENT "Running poet_init on ${BENCH_FILE}" + DEPENDS poet_init ${CMAKE_CURRENT_SOURCE_DIR}/${BENCH_FILE} + VERBATIM + COMMAND_EXPAND_LISTS + ) + + list(APPEND OUT_FILES_LIST ${OUT_FILE}) + + endforeach(BENCH_FILE ${${POET_BENCH_LIST}}) + + add_custom_target( + ${TARGET} + DEPENDS ${OUT_FILES_LIST}) + + install(FILES ${OUT_FILES_LIST} DESTINATION ${bench_install_dir}) + + # install all ADD_FILES to the same location + install(FILES ${${RT_FILES}} DESTINATION ${bench_install_dir}) +endfunction() + + +# define target name +set(BENCHTARGET benchmarks) + +add_custom_target(${BENCHTARGET} ALL) + +add_subdirectory(barite) add_subdirectory(dolo) add_subdirectory(surfex) -add_subdirectory(barite) + diff --git a/bench/barite/CMakeLists.txt b/bench/barite/CMakeLists.txt index 6c8b63aee..b51d1533b 100644 --- a/bench/barite/CMakeLists.txt +++ b/bench/barite/CMakeLists.txt @@ -1,8 +1,20 @@ -install(FILES - barite.R - barite_interp_eval.R - barite.pqi - db_barite.dat -DESTINATION - share/poet/bench/barite +# Create a list of files +set(bench_files + barite_200.R + barite_het.R ) + +set(runtime_files + barite_200_rt.R + barite_het_rt.R +) + +# add_custom_target(barite_bench DEPENDS ${bench_files} ${runtime_files}) + +ADD_BENCH_TARGET(barite_bench + bench_files + runtime_files + "barite" +) + +add_dependencies(${BENCHTARGET} barite_bench) \ No newline at end of file diff --git a/bench/barite/README.org b/bench/barite/README.org index e7b316052..86239c5de 100644 --- a/bench/barite/README.org +++ b/bench/barite/README.org @@ -18,8 +18,9 @@ mpirun -np 4 ./poet --interp barite_interp_eval.R barite_results * List of Files -- =barite.R=: POET input script for a 20x20 simulation grid -- =barite_interp_eval.R=: POET input script for a 400x200 simulation +- =barite_het.R=: POET input script with homogeneous zones for a 5x2 simulation + grid +- =barite_200.R=: POET input script for a 200x200 simulation grid - =db_barite.dat=: PHREEQC database containing the kinetic expressions for barite and celestite, stripped down from =phreeqc.dat= diff --git a/bench/barite/barite.R b/bench/barite/barite.R deleted file mode 100644 index 7d73de75b..000000000 --- a/bench/barite/barite.R +++ /dev/null @@ -1,149 +0,0 @@ -## Time-stamp: "Last modified 2024-01-12 12:39:14 delucia" -database <- normalizePath("../share/poet/bench/barite/db_barite.dat") -input_script <- normalizePath("../share/poet/bench/barite/barite.pqi") - -################################################################# -## Section 1 ## -## Grid initialization ## -################################################################# - -n <- 20 -m <- 20 - -types <- c("scratch", "phreeqc", "rds") - -init_cell <- list( - "H" = 110.0124, - "O" = 55.5087, - "Charge" = -1.216307845207e-09, - "Ba" = 1.E-12, - "Cl" = 2.E-12, - "S(6)" = 6.204727095976e-04, - "Sr" = 6.204727095976e-04, - "Barite" = 0.001, - "Celestite" = 1 -) - -grid <- list( - n_cells = c(n, m), - s_cells = c(1, 1), - type = types[1], - init_cell = as.data.frame(init_cell, check.names = FALSE), - props = names(init_cell), - database = database, - input_script = input_script -) - - -################################################################## -## Section 2 ## -## Diffusion parameters and boundary conditions ## -################################################################## - -## initial conditions - -init_diffu <- list( - "H" = 110.0124, - "O" = 55.5087, - "Charge" = -1.216307845207e-09, - "Ba" = 1.E-12, - "Cl" = 2.E-12, - "S(6)" = 6.204727095976e-04, - "Sr" = 6.204727095976e-04 -) - -injection_diff <- list( - list( - "H" = 111.0124, - "O" = 55.50622, - "Charge" = -3.336970273297e-08, - "Ba" = 0.1, - "Cl" = 0.2, - "S(6)" = 0, - "Sr" = 0) -) - -## diffusion coefficients -alpha_diffu <- c( - "H" = 1E-06, - "O" = 1E-06, - "Charge" = 1E-06, - "Ba" = 1E-06, - "Cl" = 1E-06, - "S(6)" = 1E-06, - "Sr" = 1E-06 -) - -## vecinj_inner <- list( -## l1 = c(1,20,20), -## l2 = c(2,80,80), -## l3 = c(2,60,80) -## ) - -boundary <- list( - "N" = c(1,1, rep(0, n-2)), -## "N" = rep(0, n), - "E" = rep(0, n), - "S" = rep(0, n), -# "W" = rep(0, n) - "W" = c(1,1, rep(0, n-2)) -) - -diffu_list <- names(alpha_diffu) - -vecinj <- do.call(rbind.data.frame, injection_diff) -names(vecinj) <- names(init_diffu) - -diffusion <- list( - init = as.data.frame(init_diffu, check.names = FALSE), - vecinj = vecinj, -# vecinj_inner = vecinj_inner, - vecinj_index = boundary, - alpha = alpha_diffu -) - -################################################################# -## Section 3 ## -## Chemistry module (Phreeqc) ## -################################################################# - - -## # Needed when using DHT -dht_species <- c( - "H" = 7, - "O" = 7, - "Charge" = 4, - "Ba" = 7, - "Cl" = 7, - "S(6)" = 7, - "Sr" = 7, - "Barite" = 4, - "Celestite" = 4 -) - - - -chemistry <- list( - database = database, - input_script = input_script, - dht_species = dht_species -) - -################################################################# -## Section 4 ## -## Putting all those things together ## -################################################################# - - -iterations <- 4 -dt <- 100 - -setup <- list( - grid = grid, - diffusion = diffusion, - chemistry = chemistry, - iterations = iterations, - timesteps = rep(dt, iterations), - store_result = TRUE, - out_save = seq(1, iterations) -) diff --git a/bench/barite/barite.pqi b/bench/barite/barite.pqi index 651b0e645..1fdc05e50 100644 --- a/bench/barite/barite.pqi +++ b/bench/barite/barite.pqi @@ -1,25 +1,32 @@ -SELECTED_OUTPUT - -high_precision true - -reset false - -kinetic_reactants Barite Celestite - -saturation_indices Barite Celestite SOLUTION 1 -units mol/kgw -water 1 -temperature 25 -pH 7 -pe 10.799 -Ba 0.1 -Cl 0.2 -S 1e-9 -Sr 1e-9 -KINETICS 1 -Barite --m 0.001 --parms 50. # reactive surface area --tol 1e-9 -Celestite --m 1 --parms 10.0 # reactive surface area --tol 1e-9 + units mol/kgw + water 1 + temperature 25 + pH 7 +PURE 1 + Celestite 0.0 1 END + +RUN_CELLS + -cells 1 + +COPY solution 1 2 + +KINETICS 2 + Barite + -m 0.001 + -parms 50. # reactive surface area + -tol 1e-9 + Celestite + -m 1 + -parms 10.0 # reactive surface area + -tol 1e-9 +END + +SOLUTION 3 + units mol/kgw + water 1 + temperature 25 + Ba 0.1 + Cl 0.2 +END \ No newline at end of file diff --git a/bench/barite/barite_200.R b/bench/barite/barite_200.R index 76526b00e..337f903eb 100644 --- a/bench/barite/barite_200.R +++ b/bench/barite/barite_200.R @@ -1,105 +1,39 @@ -## Time-stamp: "Last modified 2024-01-12 12:49:03 delucia" +cols <- 200 +rows <- 200 -database <- normalizePath("../share/poet/bench/barite/db_barite.dat") -input_script <- normalizePath("../share/poet/bench/barite/barite.pqi") -## database <- normalizePath("/home/work/simR/Rphree/poetsims/Sims/Hans/db_barite.dat") -## input_script <- normalizePath("/home/work/simR/Rphree/poetsims/Sims/Hans/barite.pqi") +s_cols <- 1 +s_rows <- 1 -################################################################# -## Section 1 ## -## Grid initialization ## -################################################################# +grid_def <- matrix(2, nrow = rows, ncol = cols) -n <- 200 -m <- 200 - -init_cell <- list( - "H" = 110.0124, - "O" = 55.5087, - "Charge" = -1.216307845207e-09, - "Ba" = 1.E-12, - "Cl" = 2.E-12, - "S(6)" = 6.204727095976e-04, - "Sr" = 6.204727095976e-04, - "Barite" = 0.001, - "Celestite" = 1 +# Define grid configuration for POET model +grid_setup <- list( + pqc_in_file = "./barite.pqi", + pqc_db_file = "./db_barite.dat", # Path to the database file for Phreeqc + grid_def = grid_def, # Definition of the grid, containing IDs according to the Phreeqc input script + grid_size = c(s_rows, s_cols), # Size of the grid in meters + constant_cells = c() # IDs of cells with constant concentration ) -grid <- list( - n_cells = c(n, m), - s_cells = c(1, 1), - type = "scratch", - init_cell = as.data.frame(init_cell, check.names = FALSE), - props = names(init_cell), - database = database, - input_script = input_script +bound_length <- 2 + +bound_def <- list( + "type" = rep("constant", bound_length), + "sol_id" = rep(3, bound_length), + "cell" = seq(1, bound_length) ) +homogenous_alpha <- 1e-6 -################################################################## -## Section 2 ## -## Diffusion parameters and boundary conditions ## -################################################################## - -## initial conditions - -init_diffu <- list( - "H" = 110.0124, - "O" = 55.5087, - "Charge" = -1.216307845207e-09, - "Ba" = 1.E-12, - "Cl" = 2.E-12, - "S(6)" = 6.204727095976e-04, - "Sr" = 6.204727095976e-04 +diffusion_setup <- list( + boundaries = list( + "W" = bound_def, + "N" = bound_def + ), + alpha_x = homogenous_alpha, + alpha_y = homogenous_alpha ) -injection_diff <- list( - list( - "H" = 111.0124, - "O" = 55.50622, - "Charge" = -3.336970273297e-08, - "Ba" = 0.1, - "Cl" = 0.2, - "S(6)" = 0, - "Sr" = 0) -) - -## diffusion coefficients -alpha_diffu <- c( - "H" = 1E-06, - "O" = 1E-06, - "Charge" = 1E-06, - "Ba" = 1E-06, - "Cl" = 1E-06, - "S(6)" = 1E-06, - "Sr" = 1E-06 -) - -boundary <- list( - "N" = c(1,1, rep(0, n-2)), - "E" = rep(0, n), - "S" = rep(0, n), - "W" = c(1,1, rep(0, n-2)) -) - -diffu_list <- names(alpha_diffu) - -vecinj <- do.call(rbind.data.frame, injection_diff) -names(vecinj) <- names(init_diffu) - -diffusion <- list( - init = as.data.frame(init_diffu, check.names = FALSE), - vecinj = vecinj, - vecinj_index = boundary, - alpha = alpha_diffu -) - -################################################################# -## Section 3 ## -## Chemistry module (Phreeqc) ## -################################################################# - -## DHT significant digits dht_species <- c( "H" = 7, "O" = 7, @@ -112,27 +46,13 @@ dht_species <- c( "Celestite" = 4 ) -chemistry <- list( - database = database, - input_script = input_script, - dht_species = dht_species +chemistry_setup <- list( + dht_species = dht_species ) -################################################################# -## Section 4 ## -## Putting all those things together ## -################################################################# - -iterations <- 50 -dt <- 100 - +# Define a setup list for simulation configuration setup <- list( - grid = grid, - diffusion = diffusion, - chemistry = chemistry, - iterations = iterations, - timesteps = rep(dt, iterations), - store_result = TRUE, - out_save = seq(1, iterations) - ## out_save = c(1, 5, 10, seq(50, iterations, by=50)) + Grid = grid_setup, # Parameters related to the grid structure + Diffusion = diffusion_setup, # Parameters related to the diffusion process + Chemistry = chemistry_setup ) diff --git a/bench/barite/barite_200_rt.R b/bench/barite/barite_200_rt.R new file mode 100644 index 000000000..2b4a2836d --- /dev/null +++ b/bench/barite/barite_200_rt.R @@ -0,0 +1,7 @@ +iterations <- 50 +dt <- 100 + +list( + timesteps = rep(dt, iterations), + store_result = TRUE +) \ No newline at end of file diff --git a/bench/barite/barite_het.R b/bench/barite/barite_het.R new file mode 100644 index 000000000..63807144e --- /dev/null +++ b/bench/barite/barite_het.R @@ -0,0 +1,32 @@ +grid_def <- matrix(c(2, 3), nrow = 2, ncol = 5) + +# Define grid configuration for POET model +grid_setup <- list( + pqc_in_file = "./barite_het.pqi", + pqc_db_file = "./db_barite.dat", # Path to the database file for Phreeqc + grid_def = grid_def, # Definition of the grid, containing IDs according to the Phreeqc input script + grid_size = c(ncol(grid_def), nrow(grid_def)), # Size of the grid in meters + constant_cells = c() # IDs of cells with constant concentration +) + +diffusion_setup <- list( + boundaries = list( + "W" = list( + "type" = rep("constant", nrow(grid_def)), + "sol_id" = rep(4, nrow(grid_def)), + "cell" = seq_len(nrow(grid_def)) + ) + ), + alpha_x = 1e-6, + alpha_y = matrix(runif(10, 1e-8, 1e-7), + nrow = nrow(grid_def), + ncol = ncol(grid_def) + ) +) + +# Define a setup list for simulation configuration +setup <- list( + Grid = grid_setup, # Parameters related to the grid structure + Diffusion = diffusion_setup, # Parameters related to the diffusion process + Chemistry = list() +) diff --git a/bench/barite/barite_het.pqi b/bench/barite/barite_het.pqi new file mode 100644 index 000000000..b8cf1e101 --- /dev/null +++ b/bench/barite/barite_het.pqi @@ -0,0 +1,80 @@ +## Initial: everywhere equilibrium with Celestite NB: The aqueous +## solution *resulting* from this calculation is to be used as initial +## state everywhere in the domain +SOLUTION 1 +units mol/kgw +water 1 +temperature 25 +pH 7 +pe 4 +S(6) 1e-12 +Sr 1e-12 +Ba 1e-12 +Cl 1e-12 +PURE 1 +Celestite 0.0 1 + +SAVE SOLUTION 2 # <- phreeqc keyword to store and later reuse these results +END + +RUN_CELLS + -cells 1 + +COPY solution 1 2-3 + +## Here a 5x2 domain: + + |---+---+---+---+---| + -> | 2 | 2 | 2 | 2 | 2 | + 4 |---+---+---+---+---| + -> | 3 | 3 | 3 | 3 | 3 | + |---+---+---+---+---| + +## East boundary: "injection" of solution 4. North, W, S boundaries: closed + +## Here the two distinct zones: nr 2 with kinetics Celestite (initial +## amount is 0, is then allowed to precipitate) and nr 3 with kinetic +## Celestite and Barite (both initially > 0) where the actual +## replacement takes place + +#USE SOLUTION 2 <- PHREEQC keyword to reuse the results from previous calculation +KINETICS 2 +Celestite +-m 0 # Allowed to precipitate +-parms 10.0 +-tol 1e-9 + +END + +#USE SOLUTION 2 <- PHREEQC keyword to reuse the results from previous calculation +KINETICS 3 +Barite +-m 0.001 +-parms 50. +-tol 1e-9 +Celestite +-m 1 +-parms 10.0 +-tol 1e-9 +END + +## A BaCl2 solution (nr 4) is "injected" from the left boundary: +SOLUTION 4 +units mol/kgw +pH 7 +water 1 +temp 25 +Ba 0.1 +Cl 0.2 +END +## NB: again, the *result* of the SOLUTION 4 script defines the +## concentration of all elements (+charge, tot H, tot O) + +## Ideally, in the initial state SOLUTION 1 we should not have to +## define the 4 elemental concentrations (S(6), Sr, Ba and Cl) but +## obtain them having run once the scripts with the aqueous solution +## resulting from SOLUTION 1 once with KINETICS 2 and once with +## KINETICS 3. + +RUN_CELLS + -cells 2-4 \ No newline at end of file diff --git a/bench/barite/barite_het_rt.R b/bench/barite/barite_het_rt.R new file mode 100644 index 000000000..a0b63df67 --- /dev/null +++ b/bench/barite/barite_het_rt.R @@ -0,0 +1,4 @@ +list( + timesteps = rep(50, 100), + store_result = TRUE +) \ No newline at end of file diff --git a/bench/barite/barite_interp_eval.R b/bench/barite/barite_interp_eval.R deleted file mode 100644 index 6f2b79a34..000000000 --- a/bench/barite/barite_interp_eval.R +++ /dev/null @@ -1,136 +0,0 @@ -## Time-stamp: "Last modified 2024-01-12 11:35:11 delucia" - -database <- normalizePath("../share/poet/bench/barite/db_barite.dat") -input_script <- normalizePath("../share/poet/bench/barite/barite.pqi") - -################################################################# -## Section 1 ## -## Grid initialization ## -################################################################# - -n <- 400 -m <- 200 - -types <- c("scratch", "phreeqc", "rds") - -init_cell <- list( - "H" = 110.0124, - "O" = 55.5087, - "Charge" = -1.216307845207e-09, - "Ba" = 1.E-12, - "Cl" = 2.E-12, - "S(6)" = 6.204727095976e-04, - "Sr" = 6.204727095976e-04, - "Barite" = 0.001, - "Celestite" = 1 -) - -grid <- list( - n_cells = c(n, m), - s_cells = c(1, 1), - type = types[1], - init_cell = as.data.frame(init_cell, check.names = FALSE), - props = names(init_cell), - database = database, - input_script = input_script -) - - -################################################################## -## Section 2 ## -## Diffusion parameters and boundary conditions ## -################################################################## - -## initial conditions - -init_diffu <- list( - "H" = 110.0124, - "O" = 55.5087, - "Charge" = -1.216307845207e-09, - "Ba" = 1.E-12, - "Cl" = 2.E-12, - "S(6)" = 6.204727095976e-04, - "Sr" = 6.204727095976e-04 -) - -injection_diff <- list( - list( - "H" = 111.0124, - "O" = 55.50622, - "Charge" = -3.336970273297e-08, - "Ba" = 0.1, - "Cl" = 0.2, - "S(6)" = 0, - "Sr" = 0) -) - -## diffusion coefficients -alpha_diffu <- c( - "H" = 1E-06, - "O" = 1E-06, - "Charge" = 1E-06, - "Ba" = 1E-06, - "Cl" = 1E-06, - "S(6)" = 1E-06, - "Sr" = 1E-06 -) - -boundary <- list( - "N" = c(1,1, rep(0, n-2)), - "E" = rep(0, n), - "S" = rep(0, n), - "W" = c(1,1, rep(0, n-2)) -) - -diffu_list <- names(alpha_diffu) - -vecinj <- do.call(rbind.data.frame, injection_diff) -names(vecinj) <- names(init_diffu) - -diffusion <- list( - init = as.data.frame(init_diffu, check.names = FALSE), - vecinj = vecinj, - vecinj_index = boundary, - alpha = alpha_diffu -) - -################################################################# -## Section 3 ## -## Chemistry module (Phreeqc) ## -################################################################# - -## # Needed when using DHT -dht_species <- c( - "H" = 7, - "O" = 7, - "Charge" = 4, - "Ba" = 7, - "Cl" = 7, - "S(6)" = 7, - "Sr" = 7, - "Barite" = 4, - "Celestite" = 4 -) - -chemistry <- list( - database = database, - input_script = input_script, - dht_species = dht_species -) - -################################################################# -## Section 4 ## -## Putting all those things together ## -################################################################# -iterations <- 200 -dt <- 250 - -setup <- list( - grid = grid, - diffusion = diffusion, - chemistry = chemistry, - iterations = iterations, - timesteps = rep(dt, iterations), - store_result = TRUE, - out_save = seq(1, iterations) -) diff --git a/bench/dolo/CMakeLists.txt b/bench/dolo/CMakeLists.txt index d2c0acd1e..b32d79432 100644 --- a/bench/dolo/CMakeLists.txt +++ b/bench/dolo/CMakeLists.txt @@ -1,9 +1,18 @@ -install(FILES - dolo_diffu_inner.R - dolo_diffu_inner_large.R - dolo_inner.pqi - dolo_interp_long.R - phreeqc_kin.dat -DESTINATION - share/poet/bench/dolo +set(bench_files + dolo_inner_large.R + dolo_interp.R ) + +set(runtime_files + dolo_inner_large_rt.R + dolo_interp_rt.R +) + +ADD_BENCH_TARGET( + dolo_bench + bench_files + runtime_files + "dolo" +) + +add_dependencies(${BENCHTARGET} dolo_bench) \ No newline at end of file diff --git a/bench/dolo/Eval.R b/bench/dolo/Eval.R deleted file mode 100644 index 4ac96ef1a..000000000 --- a/bench/dolo/Eval.R +++ /dev/null @@ -1,51 +0,0 @@ -## Time-stamp: "Last modified 2022-12-16 20:26:03 delucia" - -source("../../../util/data_evaluation/RFun_Eval.R") - -sd <- ReadRTSims("naaice_2d") - -sd <- ReadRTSims("Sim2D") - - -sd <- ReadRTSims("inner") - -tim <- readRDS("inner/timings.rds") - - -simtimes <- sapply(sd, "[","simtime") - -## workhorse function to be used with package "animation" -PlotAn <- function(tot, prop, grid, breaks) { - for (step in seq(1, length(tot))) { - snap <- tot[[step]]$C - time <- tot[[step]]$simtime/3600/24 - ind <- match(prop, colnames(snap)) - Plot2DCellData(snap[,ind], grid=grid, contour=FALSE, breaks=breaks, nlevels=length(breaks), scale=TRUE, main=paste0(prop," after ", time, "days")) - } -} - - -options(width=110) -library(viridis) - -Plot2DCellData(sd$iter_050$C$Cl, nx=1/100, ny=1/100, contour = TRUE, - nlevels = 12, palette = "heat.colors", - rev.palette = TRUE, scale = TRUE, main="Cl") - -Plot2DCellData(sd$iter_050$C$Dolomite, nx=100, ny=100, contour = FALSE, - nlevels = 12, palette = "heat.colors", - rev.palette = TRUE, scale = TRUE, ) - -cairo_pdf("naaice_inner_Dolo.pdf", width=8, height = 6, family="serif") -Plot2DCellData(sd$iter_100$C$Dolomite, nx=100, ny=100, contour = FALSE, - nlevels = 12, palette = "viridis", - rev.palette = TRUE, scale = TRUE, plot.axes = FALSE, - main="2D Diffusion - Dolomite after 2E+4 s (100 iterations)") -dev.off() - -cairo_pdf("naaice_inner_Mg.pdf", width=8, height = 6, family="serif") -Plot2DCellData(sd$iter_100$C$Mg, nx=100, ny=100, contour = FALSE, - nlevels = 12, palette = "terrain.colors", - rev.palette = TRUE, scale = TRUE, plot.axes=FALSE, - main="2D Diffusion - Mg after 2E+4 s (100 iterations)") -dev.off() diff --git a/bench/dolo/README.org b/bench/dolo/README.org index 67caf6e12..102b054fb 100644 --- a/bench/dolo/README.org +++ b/bench/dolo/README.org @@ -18,16 +18,13 @@ mpirun -np 4 ./poet --dht --interp dolo_interp_long.R dolo_interp_long_res * List of Files -- =dolo_diffu_inner.R=: POET input script for a 100x100 simulation - grid -- =dolo_interp_long.R=: POET input script for a 400x200 simulation +- =dolo_interp.R=: POET input script for a 400x200 simulation grid +- =dolo_diffu_inner_large.R=: POET input script for a 400x200 + simulation grid - =phreeqc_kin.dat=: PHREEQC database containing the kinetic expressions for dolomite and celestite, stripped down from =phreeqc.dat= - =dol.pqi=: PHREEQC input script for the chemical system -# - =dolo.R=: POET input script for a 20x20 simulation grid -# - =dolo_diffu_inner_large.R=: POET input script for a 400x200 -# simulation grid * Chemical system diff --git a/bench/dolo/dol.pqi b/bench/dolo/dol.pqi index 408339482..30a6674d7 100644 --- a/bench/dolo/dol.pqi +++ b/bench/dolo/dol.pqi @@ -1,35 +1,43 @@ -SELECTED_OUTPUT - -high_precision true - -reset false - -time - -soln - -temperature true - -water true - -pH - -pe - -totals C Ca Cl Mg - -kinetic_reactants Calcite Dolomite - -equilibrium O2g - SOLUTION 1 - units mol/kgw - temp 25.0 - water 1 - pH 9.91 charge - pe 4.0 - C 1.2279E-04 - Ca 1.2279E-04 - Cl 1E-12 - Mg 1E-12 + units mol/kgw + water 1 + temperature 25 + pH 7 + pe 4 PURE 1 - O2g -0.1675 10 -KINETICS 1 - Calcite - -m 0.000207 - -parms 0.0032 - -tol 1e-10 - Dolomite - -m 0.0 - -parms 0.00032 - -tol 1e-10 + Calcite 0.0 1 +END + +RUN_CELLS + -cells 1 + +COPY solution 1 2 + +PURE 2 + O2g -0.1675 10 +KINETICS 2 + Calcite + -m 0.000207 + -parms 0.05 + -tol 1e-10 + Dolomite + -m 0.0 + -parms 0.005 + -tol 1e-10 +END + +SOLUTION 3 + units mol/kgw + water 1 + temp 25 + Mg 0.001 + Cl 0.002 +END + +SOLUTION 4 + units mol/kgw + water 1 + temp 25 + Mg 0.002 + Cl 0.004 END diff --git a/bench/dolo/dolo_diffu_inner.R b/bench/dolo/dolo_diffu_inner.R deleted file mode 100644 index a16f8a49c..000000000 --- a/bench/dolo/dolo_diffu_inner.R +++ /dev/null @@ -1,190 +0,0 @@ -## Time-stamp: "Last modified 2023-08-16 17:04:42 mluebke" - -database <- normalizePath("../share/poet/bench/dolo/phreeqc_kin.dat") -input_script <- normalizePath("../share/poet/bench/dolo/dolo_inner.pqi") - -################################################################# -## Section 1 ## -## Grid initialization ## -################################################################# - -n <- 100 -m <- 100 - -types <- c("scratch", "phreeqc", "rds") - -init_cell <- list( - "H" = 110.683, - "O" = 55.3413, - "Charge" = -5.0822e-19, - "C(4)" = 1.2279E-4, - "Ca" = 1.2279E-4, - "Cl" = 0, - "Mg" = 0, - "O2g" = 0.499957, - "Calcite" = 2.07e-4, - "Dolomite" = 0 -) - -grid <- list( - n_cells = c(n, m), - s_cells = c(1, 1), - type = types[1] -) - - -################################################################## -## Section 2 ## -## Diffusion parameters and boundary conditions ## -################################################################## - -## initial conditions -init_diffu <- list( - "H" = 110.683, - "O" = 55.3413, - "Charge" = -5.0822e-19, - "C(4)" = 1.2279E-4, - "Ca" = 1.2279E-4, - "Cl" = 0, - "Mg" = 0 -) - -## diffusion coefficients -alpha_diffu <- c( - "H" = 1E-6, - "O" = 1E-6, - "Charge" = 1E-6, - "C(4)" = 1E-6, - "Ca" = 1E-6, - "Cl" = 1E-6, - "Mg" = 1E-6 -) - -## list of boundary conditions/inner nodes -vecinj_diffu <- list( - list( - "H" = 110.683, - "O" = 55.3413, - "Charge" = 1.90431e-16, - "C(4)" = 0, - "Ca" = 0, - "Cl" = 0.002, - "Mg" = 0.001 - ), - list( - "H" = 110.683, - "O" = 55.3413, - "Charge" = 1.90431e-16, - "C(4)" = 0, - "Ca" = 0.0, - "Cl" = 0.004, - "Mg" = 0.002 - ) -) - -vecinj_inner <- list( - l1 = c(1, 20, 20), - l2 = c(2, 80, 80), - l3 = c(2, 60, 80) -) - -boundary <- list( - # "N" = c(1, rep(0, n-1)), - "N" = rep(0, n), - "E" = rep(0, n), - "S" = rep(0, n), - "W" = rep(0, n) -) - -diffu_list <- names(alpha_diffu) - -vecinj <- do.call(rbind.data.frame, vecinj_diffu) -names(vecinj) <- names(init_diffu) - -diffusion <- list( - init = as.data.frame(init_diffu, check.names = FALSE), - vecinj = vecinj, - vecinj_inner = vecinj_inner, - vecinj_index = boundary, - alpha = alpha_diffu -) - -################################################################# -## Section 3 ## -## Chemistry module (Phreeqc) ## -################################################################# - - -## # Needed when using DHT -dht_species <- c( - "H" = 10, - "O" = 10, - "Charge" = 3, - "C(4)" = 5, - "Ca" = 5, - "Cl" = 5, - "Mg" = 5, - "Calcite" = 5, - "Dolomite" = 5 -) - -check_sign_cal_dol_dht <- function(old, new) { - if ((old["Calcite"] == 0) != (new["Calcite"] == 0)) { - return(TRUE) - } - if ((old["Dolomite"] == 0) != (new["Dolomite"] == 0)) { - return(TRUE) - } - return(FALSE) -} - -fuzz_input_dht_keys <- function(input) { - return(input[names(dht_species)]) -} - -check_sign_cal_dol_interp <- function(to_interp, data_set) { - data_set <- as.data.frame(do.call(rbind, data_set), check.names = FALSE, optional = TRUE) - names(data_set) <- names(dht_species) - cal <- (data_set$Calcite == 0) == (to_interp["Calcite"] == 0) - dol <- (data_set$Dolomite == 0) == (to_interp["Dolomite"] == 0) - - cal_dol_same_sig <- cal == dol - return(rev(which(!cal_dol_same_sig))) -} - -check_neg_cal_dol <- function(result) { - neg_sign <- (result["Calcite"] <- 0) || (result["Dolomite"] < 0) - return(any(neg_sign)) -} - -hooks <- list( - dht_fill = check_sign_cal_dol_dht, - dht_fuzz = fuzz_input_dht_keys, - interp_pre_func = check_sign_cal_dol_interp, - interp_post_func = check_neg_cal_dol -) - -chemistry <- list( - database = database, - input_script = input_script, - dht_species = dht_species, - hooks = hooks -) - -################################################################# -## Section 4 ## -## Putting all those things together ## -################################################################# - - -iterations <- 10 -dt <- 200 - -setup <- list( - grid = grid, - diffusion = diffusion, - chemistry = chemistry, - iterations = iterations, - timesteps = rep(dt, iterations), - store_result = TRUE -) diff --git a/bench/dolo/dolo_diffu_inner_large.R b/bench/dolo/dolo_diffu_inner_large.R deleted file mode 100644 index c9c70c6a4..000000000 --- a/bench/dolo/dolo_diffu_inner_large.R +++ /dev/null @@ -1,190 +0,0 @@ -## Time-stamp: "Last modified 2023-08-16 17:05:04 mluebke" - -database <- normalizePath("../share/poet/bench/dolo/phreeqc_kin.dat") -input_script <- normalizePath("../share/poet/bench/dolo/dolo_inner.pqi") - -################################################################# -## Section 1 ## -## Grid initialization ## -################################################################# - -n <- 2000 -m <- 1000 - -types <- c("scratch", "phreeqc", "rds") - -init_cell <- list( - "H" = 110.683, - "O" = 55.3413, - "Charge" = -5.0822e-19, - "C" = 1.2279E-4, - "Ca" = 1.2279E-4, - "Cl" = 0, - "Mg" = 0, - "O2g" = 0.499957, - "Calcite" = 2.07e-4, - "Dolomite" = 0 -) - -grid <- list( - n_cells = c(n, m), - s_cells = c(20, 10), - type = types[1] -) - - -################################################################## -## Section 2 ## -## Diffusion parameters and boundary conditions ## -################################################################## - -## initial conditions -init_diffu <- list( - "H" = 110.683, - "O" = 55.3413, - "Charge" = -5.0822e-19, - "C(4)" = 1.2279E-4, - "Ca" = 1.2279E-4, - "Cl" = 0, - "Mg" = 0 -) - -## diffusion coefficients -alpha_diffu <- c( - "H" = 1E-6, - "O" = 1E-6, - "Charge" = 1E-6, - "C(4)" = 1E-6, - "Ca" = 1E-6, - "Cl" = 1E-6, - "Mg" = 1E-6 -) - -## list of boundary conditions/inner nodes -vecinj_diffu <- list( - list( - "H" = 110.683, - "O" = 55.3413, - "Charge" = 1.90431e-16, - "C(4)" = 0, - "Ca" = 0, - "Cl" = 0.002, - "Mg" = 0.001 - ), - list( - "H" = 110.683, - "O" = 55.3413, - "Charge" = 1.90431e-16, - "C(4)" = 0, - "Ca" = 0.0, - "Cl" = 0.004, - "Mg" = 0.002 - ) -) - -vecinj_inner <- list( - l1 = c(1, 400, 200), - l2 = c(2, 1400, 800), - l3 = c(2, 1600, 800) -) - -boundary <- list( - # "N" = c(1, rep(0, n-1)), - "N" = rep(0, n), - "E" = rep(0, m), - "S" = rep(0, n), - "W" = rep(0, m) -) - -diffu_list <- names(alpha_diffu) - -vecinj <- do.call(rbind.data.frame, vecinj_diffu) -names(vecinj) <- names(init_diffu) - -diffusion <- list( - init = as.data.frame(init_diffu, check.names = FALSE), - vecinj = vecinj, - vecinj_inner = vecinj_inner, - vecinj_index = boundary, - alpha = alpha_diffu -) - -################################################################# -## Section 3 ## -## Chemistry module (Phreeqc) ## -################################################################# - -## # Needed when using DHT -dht_species <- c( - "H" = 10, - "O" = 10, - "Charge" = 3, - "C(4)" = 5, - "Ca" = 5, - "Cl" = 5, - "Mg" = 5, - "Calcite" = 5, - "Dolomite" = 5 -) - -check_sign_cal_dol_dht <- function(old, new) { - if ((old["Calcite"] == 0) != (new["Calcite"] == 0)) { - return(TRUE) - } - if ((old["Dolomite"] == 0) != (new["Dolomite"] == 0)) { - return(TRUE) - } - return(FALSE) -} - -fuzz_input_dht_keys <- function(input) { - return(input[names(dht_species)]) -} - -check_sign_cal_dol_interp <- function(to_interp, data_set) { - data_set <- as.data.frame(do.call(rbind, data_set), check.names = FALSE, optional = TRUE) - names(data_set) <- names(dht_species) - cal <- (data_set$Calcite == 0) == (to_interp["Calcite"] == 0) - dol <- (data_set$Dolomite == 0) == (to_interp["Dolomite"] == 0) - - cal_dol_same_sig <- cal == dol - return(rev(which(!cal_dol_same_sig))) -} - -check_neg_cal_dol <- function(result) { - neg_sign <- (result["Calcite"] <- 0) || (result["Dolomite"] < 0) - return(any(neg_sign)) -} - -hooks <- list( - dht_fill = check_sign_cal_dol_dht, - dht_fuzz = fuzz_input_dht_keys, - interp_pre_func = check_sign_cal_dol_interp, - interp_post_func = check_neg_cal_dol -) - -chemistry <- list( - database = database, - input_script = input_script, - dht_species = dht_species, - hooks = hooks -) - -################################################################# -## Section 4 ## -## Putting all those things together ## -################################################################# - - -iterations <- 500 -dt <- 50 - -setup <- list( - grid = grid, - diffusion = diffusion, - chemistry = chemistry, - iterations = iterations, - timesteps = rep(dt, iterations), - store_result = TRUE, - out_save = seq(5, iterations, by = 5) -) diff --git a/bench/dolo/dolo_inner.pqi b/bench/dolo/dolo_inner.pqi deleted file mode 100644 index 99e3a2f89..000000000 --- a/bench/dolo/dolo_inner.pqi +++ /dev/null @@ -1,28 +0,0 @@ -SELECTED_OUTPUT - -high_precision true - -reset false - -kinetic_reactants Calcite Dolomite - -equilibrium O2g - -SOLUTION 1 - units mol/kgw - temp 25.0 - water 1 - pH 9.91 charge - pe 4.0 - C 1.2279E-04 - Ca 1.2279E-04 - Cl 1E-12 - Mg 1E-12 -PURE 1 - O2g -0.1675 10 -KINETICS 1 - Calcite - -m 0.00020 - -parms 0.05 - -tol 1e-10 - Dolomite - -m 0.0 - -parms 0.005 - -tol 1e-10 -END diff --git a/bench/dolo/dolo_inner.rds b/bench/dolo/dolo_inner.rds new file mode 100644 index 000000000..18740ebfc Binary files /dev/null and b/bench/dolo/dolo_inner.rds differ diff --git a/bench/dolo/dolo_inner_large.R b/bench/dolo/dolo_inner_large.R new file mode 100644 index 000000000..155f75c85 --- /dev/null +++ b/bench/dolo/dolo_inner_large.R @@ -0,0 +1,115 @@ +rows <- 2000 +cols <- 1000 + +grid_def <- matrix(2, nrow = rows, ncol = cols) + +# Define grid configuration for POET model +grid_setup <- list( + pqc_in_file = "./dol.pqi", + pqc_db_file = "./phreeqc_kin.dat", # Path to the database file for Phreeqc + grid_def = grid_def, # Definition of the grid, containing IDs according to the Phreeqc input script + grid_size = c(cols, rows) / 100, # Size of the grid in meters + constant_cells = c() # IDs of cells with constant concentration +) + +bound_size <- 2 + +diffusion_setup <- list( + inner_boundaries = list( + "row" = c(400, 1400, 1600), + "col" = c(200, 800, 800), + "sol_id" = c(3, 4, 4) + ), + alpha_x = 1e-6, + alpha_y = 1e-6 +) + +check_sign_cal_dol_dht <- function(old, new) { + if ((old["Calcite"] == 0) != (new["Calcite"] == 0)) { + return(TRUE) + } + if ((old["Dolomite"] == 0) != (new["Dolomite"] == 0)) { + return(TRUE) + } + return(FALSE) +} + +fuzz_input_dht_keys <- function(input) { + dht_species <- c( + "H" = 3, + "O" = 3, + "Charge" = 3, + "C(4)" = 6, + "Ca" = 6, + "Cl" = 3, + "Mg" = 5, + "Calcite" = 4, + "Dolomite" = 4 + ) + return(input[names(dht_species)]) +} + +check_sign_cal_dol_interp <- function(to_interp, data_set) { + dht_species <- c( + "H" = 3, + "O" = 3, + "Charge" = 3, + "C(4)" = 6, + "Ca" = 6, + "Cl" = 3, + "Mg" = 5, + "Calcite" = 4, + "Dolomite" = 4 + ) + data_set <- as.data.frame(do.call(rbind, data_set), check.names = FALSE, optional = TRUE) + names(data_set) <- names(dht_species) + cal <- (data_set$Calcite == 0) == (to_interp["Calcite"] == 0) + dol <- (data_set$Dolomite == 0) == (to_interp["Dolomite"] == 0) + + cal_dol_same_sig <- cal == dol + return(rev(which(!cal_dol_same_sig))) +} + +check_neg_cal_dol <- function(result) { + neg_sign <- (result["Calcite"] < 0) || (result["Dolomite"] < 0) + return(neg_sign) +} + +# Optional when using Interpolation (example with less key species and custom +# significant digits) + +pht_species <- c( + "C(4)" = 3, + "Ca" = 3, + "Mg" = 2, + "Calcite" = 2, + "Dolomite" = 2 +) + +chemistry_setup <- list( + dht_species = c( + "H" = 3, + "O" = 3, + "Charge" = 3, + "C(4)" = 6, + "Ca" = 6, + "Cl" = 3, + "Mg" = 5, + "Calcite" = 4, + "Dolomite" = 4 + ), + pht_species = pht_species, + hooks = list( + dht_fill = check_sign_cal_dol_dht, + dht_fuzz = fuzz_input_dht_keys, + interp_pre = check_sign_cal_dol_interp, + interp_post = check_neg_cal_dol + ) +) + +# Define a setup list for simulation configuration +setup <- list( + Grid = grid_setup, # Parameters related to the grid structure + Diffusion = diffusion_setup, # Parameters related to the diffusion process + Chemistry = chemistry_setup # Parameters related to the chemistry process +) diff --git a/bench/dolo/dolo_inner_large_rt.R b/bench/dolo/dolo_inner_large_rt.R new file mode 100644 index 000000000..27dda1f45 --- /dev/null +++ b/bench/dolo/dolo_inner_large_rt.R @@ -0,0 +1,10 @@ +iterations <- 500 +dt <- 50 + +out_save <- seq(5, iterations, by = 5) + +list( + timesteps = rep(dt, iterations), + store_result = TRUE, + out_save = out_save +) diff --git a/bench/dolo/dolo_interp.R b/bench/dolo/dolo_interp.R new file mode 100644 index 000000000..f48527d00 --- /dev/null +++ b/bench/dolo/dolo_interp.R @@ -0,0 +1,131 @@ +rows <- 400 +cols <- 200 + +grid_def <- matrix(2, nrow = rows, ncol = cols) + +# Define grid configuration for POET model +grid_setup <- list( + pqc_in_file = "./dol.pqi", + pqc_db_file = "./phreeqc_kin.dat", # Path to the database file for Phreeqc + grid_def = grid_def, # Definition of the grid, containing IDs according to the Phreeqc input script + grid_size = c(2.5, 5), # Size of the grid in meters + constant_cells = c() # IDs of cells with constant concentration +) + +bound_def_we <- list( + "type" = rep("constant", rows), + "sol_id" = rep(1, rows), + "cell" = seq(1, rows) +) + +bound_def_ns <- list( + "type" = rep("constant", cols), + "sol_id" = rep(1, cols), + "cell" = seq(1, cols) +) + +diffusion_setup <- list( + boundaries = list( + "W" = bound_def_we, + "E" = bound_def_we, + "N" = bound_def_ns, + "S" = bound_def_ns + ), + inner_boundaries = list( + "row" = floor(rows / 2), + "col" = floor(cols / 2), + "sol_id" = c(3) + ), + alpha_x = 1e-6, + alpha_y = 1e-6 +) + +check_sign_cal_dol_dht <- function(old, new) { + if ((old["Calcite"] == 0) != (new["Calcite"] == 0)) { + return(TRUE) + } + if ((old["Dolomite"] == 0) != (new["Dolomite"] == 0)) { + return(TRUE) + } + return(FALSE) +} + +fuzz_input_dht_keys <- function(input) { + dht_species <- c( + "H" = 3, + "O" = 3, + "Charge" = 3, + "C(4)" = 6, + "Ca" = 6, + "Cl" = 3, + "Mg" = 5, + "Calcite" = 4, + "Dolomite" = 4 + ) + return(input[names(dht_species)]) +} + +check_sign_cal_dol_interp <- function(to_interp, data_set) { + dht_species <- c( + "H" = 3, + "O" = 3, + "Charge" = 3, + "C(4)" = 6, + "Ca" = 6, + "Cl" = 3, + "Mg" = 5, + "Calcite" = 4, + "Dolomite" = 4 + ) + data_set <- as.data.frame(do.call(rbind, data_set), check.names = FALSE, optional = TRUE) + names(data_set) <- names(dht_species) + cal <- (data_set$Calcite == 0) == (to_interp["Calcite"] == 0) + dol <- (data_set$Dolomite == 0) == (to_interp["Dolomite"] == 0) + + cal_dol_same_sig <- cal == dol + return(rev(which(!cal_dol_same_sig))) +} + +check_neg_cal_dol <- function(result) { + neg_sign <- (result["Calcite"] < 0) || (result["Dolomite"] < 0) + return(neg_sign) +} + +# Optional when using Interpolation (example with less key species and custom +# significant digits) + +pht_species <- c( + "C(4)" = 3, + "Ca" = 3, + "Mg" = 2, + "Calcite" = 2, + "Dolomite" = 2 +) + +chemistry_setup <- list( + dht_species = c( + "H" = 3, + "O" = 3, + "Charge" = 3, + "C(4)" = 6, + "Ca" = 6, + "Cl" = 3, + "Mg" = 5, + "Calcite" = 4, + "Dolomite" = 4 + ), + pht_species = pht_species, + hooks = list( + dht_fill = check_sign_cal_dol_dht, + dht_fuzz = fuzz_input_dht_keys, + interp_pre = check_sign_cal_dol_interp, + interp_post = check_neg_cal_dol + ) +) + +# Define a setup list for simulation configuration +setup <- list( + Grid = grid_setup, # Parameters related to the grid structure + Diffusion = diffusion_setup, # Parameters related to the diffusion process + Chemistry = chemistry_setup # Parameters related to the chemistry process +) diff --git a/bench/dolo/dolo_interp_long.R b/bench/dolo/dolo_interp_long.R deleted file mode 100644 index 7acca4298..000000000 --- a/bench/dolo/dolo_interp_long.R +++ /dev/null @@ -1,204 +0,0 @@ -## Time-stamp: "Last modified 2023-08-16 14:57:25 mluebke" - -database <- normalizePath("../share/poet/bench/dolo/phreeqc_kin.dat") -input_script <- normalizePath("../share/poet/bench/dolo/dolo_inner.pqi") - -################################################################# -## Section 1 ## -## Grid initialization ## -################################################################# - -n <- 400 -m <- 200 - -types <- c("scratch", "phreeqc", "rds") - -init_cell <- list( - "H" = 110.683, - "O" = 55.3413, - "Charge" = -5.0822e-19, - "C" = 1.2279E-4, - "Ca" = 1.2279E-4, - "Cl" = 0, - "Mg" = 0, - "O2g" = 0.499957, - "Calcite" = 2.07e-4, - "Dolomite" = 0 -) - -grid <- list( - n_cells = c(n, m), - s_cells = c(5, 2.5), - type = types[1] -) - - -################################################################## -## Section 2 ## -## Diffusion parameters and boundary conditions ## -################################################################## - -## initial conditions -init_diffu <- list( - "H" = 1.110124E+02, - "O" = 5.550833E+01, - "Charge" = -1.216307659761E-09, - "C(4)" = 1.230067028174E-04, - "Ca" = 1.230067028174E-04, - "Cl" = 0, - "Mg" = 0 -) - -## diffusion coefficients -alpha_diffu <- c( - "H" = 1E-6, - "O" = 1E-6, - "Charge" = 1E-6, - "C(4)" = 1E-6, - "Ca" = 1E-6, - "Cl" = 1E-6, - "Mg" = 1E-6 -) - -## list of boundary conditions/inner nodes -vecinj_diffu <- list( - list( - "H" = 1.110124E+02, - "O" = 5.550796E+01, - "Charge" = -3.230390327801E-08, - "C(4)" = 0, - "Ca" = 0, - "Cl" = 0.002, - "Mg" = 0.001 - ), - list( - "H" = 110.683, - "O" = 55.3413, - "Charge" = 1.90431e-16, - "C(4)" = 0, - "Ca" = 0.0, - "Cl" = 0.004, - "Mg" = 0.002 - ), - init_diffu -) - -vecinj_inner <- list( - l1 = c(1, floor(n / 2), floor(m / 2)) - # l2 = c(2,1400,800), - # l3 = c(2,1600,800) -) - -boundary <- list( - # "N" = c(1, rep(0, n-1)), - "N" = rep(3, n), - "E" = rep(3, m), - "S" = rep(3, n), - "W" = rep(3, m) -) - -diffu_list <- names(alpha_diffu) - -vecinj <- do.call(rbind.data.frame, vecinj_diffu) -names(vecinj) <- names(init_diffu) - -diffusion <- list( - init = as.data.frame(init_diffu, check.names = FALSE), - vecinj = vecinj, - vecinj_inner = vecinj_inner, - vecinj_index = boundary, - alpha = alpha_diffu -) - -################################################################# -## Section 3 ## -## Chemistry module (Phreeqc) ## -################################################################# - - -## # optional when using DHT -dht_species <- c( - "H" = 3, - "O" = 3, - "Charge" = 3, - "C(4)" = 6, - "Ca" = 6, - "Cl" = 3, - "Mg" = 5, - "Calcite" = 4, - "Dolomite" = 4 -) - -## # Optional when using Interpolation (example with less key species and custom -## # significant digits) - -# pht_species <- c( -# "C(4)" = 3, -# "Ca" = 3, -# "Mg" = 2, -# "Calcite" = 2, -# "Dolomite" = 2 -# ) - -check_sign_cal_dol_dht <- function(old, new) { - if ((old["Calcite"] == 0) != (new["Calcite"] == 0)) { - return(TRUE) - } - if ((old["Dolomite"] == 0) != (new["Dolomite"] == 0)) { - return(TRUE) - } - return(FALSE) -} - -fuzz_input_dht_keys <- function(input) { - return(input[names(dht_species)]) -} - -check_sign_cal_dol_interp <- function(to_interp, data_set) { - data_set <- as.data.frame(do.call(rbind, data_set), check.names = FALSE, optional = TRUE) - names(data_set) <- names(dht_species) - cal <- (data_set$Calcite == 0) == (to_interp["Calcite"] == 0) - dol <- (data_set$Dolomite == 0) == (to_interp["Dolomite"] == 0) - - cal_dol_same_sig <- cal == dol - return(rev(which(!cal_dol_same_sig))) -} - -check_neg_cal_dol <- function(result) { - neg_sign <- (result["Calcite"] < 0) || (result["Dolomite"] < 0) - return(neg_sign) -} - -hooks <- list( - dht_fill = check_sign_cal_dol_dht, - dht_fuzz = fuzz_input_dht_keys, - interp_pre_func = check_sign_cal_dol_interp, - interp_post_func = check_neg_cal_dol -) - -chemistry <- list( - database = database, - input_script = input_script, - dht_species = dht_species, - hooks = hooks - # pht_species = pht_species -) - -################################################################# -## Section 4 ## -## Putting all those things together ## -################################################################# - - -iterations <- 20000 -dt <- 200 - -setup <- list( - grid = grid, - diffusion = diffusion, - chemistry = chemistry, - iterations = iterations, - timesteps = rep(dt, iterations), - store_result = TRUE, - out_save = c(1, seq(50, iterations, by = 50)) -) diff --git a/bench/dolo/dolo_interp_rt.R b/bench/dolo/dolo_interp_rt.R new file mode 100644 index 000000000..417cd5b8d --- /dev/null +++ b/bench/dolo/dolo_interp_rt.R @@ -0,0 +1,10 @@ +iterations <- 20000 +dt <- 200 + +out_save <- seq(50, iterations, by = 50) + +list( + timesteps = rep(dt, iterations), + store_result = TRUE, + out_save = out_save +) diff --git a/bench/surfex/CMakeLists.txt b/bench/surfex/CMakeLists.txt index c5a605c6c..baf23bcb0 100644 --- a/bench/surfex/CMakeLists.txt +++ b/bench/surfex/CMakeLists.txt @@ -1,9 +1,20 @@ -install(FILES - ExBase.pqi - ex.R - surfex.R - SurfExBase.pqi - SMILE_2021_11_01_TH.dat -DESTINATION - share/poet/bench/surfex +set(bench_files + # surfex.R + # ex.R + PoetEGU_surfex_500.R ) + +set(runtime_files + # surfex_rt.R + # ex_rt.R + PoetEGU_surfex_500_rt.R +) + +ADD_BENCH_TARGET( + surfex_bench + bench_files + runtime_files + "surfex" +) + +add_dependencies(${BENCHTARGET} surfex_bench) \ No newline at end of file diff --git a/bench/surfex/ExBase.pqi b/bench/surfex/ExBase.pqi index af9c945a3..908488f1c 100644 --- a/bench/surfex/ExBase.pqi +++ b/bench/surfex/ExBase.pqi @@ -37,3 +37,27 @@ EXCHANGE 1 Z 0.0012585 Y 0.0009418 END + +SOLUTION 2 +temp 13 +units mol/kgw + +C(-4) 2.92438561098248e-21 +C(4) 2.65160558871092e-06 +Ca 2.89001071336443e-05 +Cl 0.000429291158114428 +Fe(2) 1.90823391198114e-07 +Fe(3) 3.10832423034763e-12 +H(0) 2.7888235127385e-15 +K 2.5301787e-06 +Mg 2.31391999937907e-05 +Na 0.00036746969 +S(-2) 1.01376078438546e-14 +S(2) 1.42247026981542e-19 +S(4) 9.49422092568557e-18 +S(6) 2.19812504654191e-05 +Sr 6.01218519999999e-07 +U(4) 4.82255946569383e-12 +U(5) 5.49050615347901e-13 +U(6) 1.32462838991902e-09 +END diff --git a/bench/surfex/PoetEGU_surfex_500.R b/bench/surfex/PoetEGU_surfex_500.R new file mode 100644 index 000000000..d6b32fd78 --- /dev/null +++ b/bench/surfex/PoetEGU_surfex_500.R @@ -0,0 +1,40 @@ +rows <- 500 +cols <- 200 + +grid_left <- matrix(1, nrow = rows, ncol = cols/2) +grid_rght <- matrix(2, nrow = rows, ncol = cols/2) +grid_def <- cbind(grid_left, grid_rght) + + +# Define grid configuration for POET model +grid_setup <- list( + pqc_in_file = "./SurfexEGU.pqi", + pqc_db_file = "./SMILE_2021_11_01_TH.dat", # Path to the database file for Phreeqc + grid_def = grid_def, # Definition of the grid, containing IDs according to the Phreeqc input script + grid_size = c(10, 4), # Size of the grid in meters + constant_cells = c() # IDs of cells with constant concentration +) + +bound_def <- list( + "type" = rep("constant", cols), + "sol_id" = rep(3, cols), + "cell" = seq(1, cols) +) + +diffusion_setup <- list( + boundaries = list( + "N" = bound_def + ), + alpha_x = matrix(runif(rows*cols))*1e-8, + alpha_y = matrix(runif(rows*cols))*1e-9## ,1e-10 +) + + +chemistry_setup <- list() + +# Define a setup list for simulation configuration +setup <- list( + Grid = grid_setup, # Parameters related to the grid structure + Diffusion = diffusion_setup, # Parameters related to the diffusion process + Chemistry = chemistry_setup # Parameters related to the chemistry process +) diff --git a/bench/surfex/PoetEGU_surfex_500_rt.R b/bench/surfex/PoetEGU_surfex_500_rt.R new file mode 100644 index 000000000..d91e17b3a --- /dev/null +++ b/bench/surfex/PoetEGU_surfex_500_rt.R @@ -0,0 +1,11 @@ +iterations <- 200 +dt <- 1000 + +out_save <- c(1, 2, seq(5, iterations, by=5)) +## out_save <- seq(1, iterations) + +list( + timesteps = rep(dt, iterations), + store_result = TRUE, + out_save = out_save +) diff --git a/bench/surfex/README.org b/bench/surfex/README.org index b765f4ef5..c99f1184e 100644 --- a/bench/surfex/README.org +++ b/bench/surfex/README.org @@ -20,7 +20,7 @@ mpirun -np 4 ./poet surfex.R surfex_res - =ex.R=: POET input script for a 100x100 simulation grid, only exchange - =ExBase.pqi=: PHREEQC input script for the =ex.R= model -- =surfex.R=: POET input script for a 100x100 simulation grid +- =surfex.R=: POET input script for a 1000x1000 simulation grid considering both cation exchange and surface complexation - =SurfExBase.pqi=: PHREEQC input script for the =surfex.R= model - =SMILE_2021_11_01_TH.dat=: PHREEQC database containing the diff --git a/bench/surfex/SurfExBase.pqi b/bench/surfex/SurfExBase.pqi index a820d70e0..2598b72b2 100644 --- a/bench/surfex/SurfExBase.pqi +++ b/bench/surfex/SurfExBase.pqi @@ -54,3 +54,27 @@ EXCHANGE 1 Z 0.0012585 Y 0.0009418 END + +SOLUTION 2 +temp 13 +units mol/kgw + +C(-4) 2.92438561098248e-21 +C(4) 2.65160558871092e-06 +Ca 2.89001071336443e-05 +Cl 0.000429291158114428 +Fe(2) 1.90823391198114e-07 +Fe(3) 3.10832423034763e-12 +H(0) 2.7888235127385e-15 +K 2.5301787e-06 +Mg 2.31391999937907e-05 +Na 0.00036746969 +S(-2) 1.01376078438546e-14 +S(2) 1.42247026981542e-19 +S(4) 9.49422092568557e-18 +S(6) 2.19812504654191e-05 +Sr 6.01218519999999e-07 +U(4) 4.82255946569383e-12 +U(5) 5.49050615347901e-13 +U(6) 1.32462838991902e-09 +END \ No newline at end of file diff --git a/bench/surfex/SurfexEGU.pqi b/bench/surfex/SurfexEGU.pqi new file mode 100644 index 000000000..a2b3342b2 --- /dev/null +++ b/bench/surfex/SurfexEGU.pqi @@ -0,0 +1,108 @@ +## Time-stamp: "Last modified 2024-04-12 10:59:59 delucia" +## KNOBS +## -logfile false +## -iterations 10000 +## -convergence_tolerance 1E-12 +## -step_size 2 +## -pe_step_size 2 + +SOLUTION 1 ## Porewater composition Opalinus Clay, WITHOUT radionuclides, AFTER EQUI_PHASES + pe -2.627 ## Eh = -227 mV, Value from Bossart & Thury (2008)-> PC borehole measurement 2003, Eh still decreasing + density 1.01583 ## kg/dm³ = g/cm³ + temp 13 ## mean temperature Mont Terri, Bossart & Thury (2008), calculations performed for 25°C + units mol/kgw + ## Mean composition + pH 7.064 + Na 2.763e-01 + Cl 3.228e-01 charge + S(6) 1.653e-02 as SO4 + Ca 2.173e-02 + Mg 1.740e-02 + K 1.902e-03 + Sr 4.520e-04 + Fe 1.435e-04 + U 2.247e-09 + +SURFACE 1 Opalinus Clay, clay minerals + ## calculated with rho_b=2.2903 kg/dm³, poro=0.1662 + ## 1 dm³ = 13.565641 kg_sed/kg_pw + -equil 1 ## equilibrate with solution 1 + -sites_units density ## set unit for binding site density to sites/nm2 + -donnan 4.9e-10 ## calculated after Wigger & Van Loon (2018) for ionic strength after equilibration with minerales for pCO2=2.2 log10 bar + + # surface density SSA (m2/g) mass (g/kgw) + Kln_aOH 1.155 11. 3798.4 ## Kaolinite 28 wt% (aluminol and silanol sites) + Kln_siOH 1.155 + Ill_sOH 0.05 100. 4205.35 ## Illite 31 wt% (weak und strong binding sites) + Ill_wOH 2.26 ## 2 % strong binding sites + Mll_sOH 0.05 100. 813.94 ## Montmorillonite = smektite = 6 wt% (weak und strong binding sites) + Mll_wOH 2.26 ## 2 % strong binding sites + +EXCHANGE 1 Exchanger, only illite+montmorillonite + ## Illite = 0.225 eq/kg_rock, Montmorillonit = 0.87 eq/kg_rock + -equil 1 ## equilibrate with solution 1 + Z 0.9462 ## = Illite + Y 0.70813 ## = Montmorillonite +END + +SOLUTION 2 ## Porewater composition Opalinus Clay, WITHOUT radionuclides, AFTER EQUI_PHASES + pe -2.627 ## Eh = -227 mV, Value from Bossart & Thury (2008)-> PC borehole measurement 2003, Eh still decreasing + density 1.01583 ## kg/dm³ = g/cm³ + temp 13 ## mean temperature Mont Terri, Bossart & Thury (2008), calculations performed for 25°C + units mol/kgw + ## Mean composition + pH 7.064 + Na 2.763e-01 + Cl 3.228e-01 charge + S(6) 1.653e-02 as SO4 + Ca 2.173e-02 + Mg 1.740e-02 + K 1.902e-03 + Sr 4.520e-04 + Fe 1.435e-04 + U 2.247e-09 + +SURFACE 2 Opalinus Clay, clay minerals + -equil 2 ## equilibrate with solution 2 + -sites_units density ## set unit for binding site density to + ## sites/nm2 + -donnan 4.9e-10 ## calculated after Wigger & Van Loon (2018) + ## for ionic strength after equilibration + ## with minerales for pCO2=2.2 log10 bar + + ## surface density SSA (m2/g) mass (g/kgw) + Kln_aOH 1.155 11. 2798.4 ## Kaolinite 28 wt% (aluminol and silanol sites) + Kln_siOH 1.155 + Ill_sOH 0.05 100. 1205.35 ## Illite 31 wt% (weak und strong binding sites) + Ill_wOH 2.26 ## 2 % strong binding sites + Mll_sOH 0.05 100. 113.94 ## Montmorillonite = smektite = 6 wt% (weak und strong binding sites) + Mll_wOH 2.26 ## 2 % strong binding sites + +EXCHANGE 2 Exchanger, only illite+montmorillonite + ## Illite = 0.225 eq/kg_rock, Montmorillonit = 0.87 eq/kg_rock + -equil 2 ## equilibrate with solution 1 + Z 0.5 ## = Illite + Y 0.2 ## = Montmorillonite +END + +SOLUTION 3 + pe -2.627 ## Eh = -227 mV, Value from Bossart & Thury (2008)-> PC borehole measurement 2003, Eh still decreasing + density 1.01583 ## kg/dm³ = g/cm³ + temp 13 ## mean temperature Mont Terri, Bossart & Thury (2008), calculations performed for 25°C + units mol/kgw + ## Mean composition + pH 7.064 + Na 3.763e-01 + Cl 4.228e-01 charge + S(6) 1.653e-02 as SO4 + Ca 2.173e-02 + Mg 1.740e-02 + K 1.902e-03 + Sr 4.520e-04 + Fe 1.435e-04 + U 1e-6 + C 1.991e-03 +END + +RUN_CELLS +END diff --git a/bench/surfex/ex.R b/bench/surfex/ex.R index 13feea1d8..2ba5399bb 100644 --- a/bench/surfex/ex.R +++ b/bench/surfex/ex.R @@ -1,140 +1,37 @@ -## Time-stamp: "Last modified 2023-08-02 13:59:35 mluebke" +rows <- 100 +cols <- 100 -database <- normalizePath("./SMILE_2021_11_01_TH.dat") -input_script <- normalizePath("./ExBase.pqi") +grid_def <- matrix(1, nrow = rows, ncol = cols) -cat(paste(":: R This is a test 1\n")) +# Define grid configuration for POET model +grid_setup <- list( + pqc_in_file = "./SurfExBase.pqi", + pqc_db_file = "./SMILE_2021_11_01_TH.dat", # Path to the database file for Phreeqc + grid_def = grid_def, # Definition of the grid, containing IDs according to the Phreeqc input script + grid_size = c(1, 1), # Size of the grid in meters + constant_cells = c() # IDs of cells with constant concentration +) -################################################################# -## Section 1 ## -## Grid initialization ## -################################################################# +bound_def <- list( + "type" = rep("constant", cols), + "sol_id" = rep(2, cols), + "cell" = seq(1, cols) +) -n <- 100 -m <- 100 - -types <- c("scratch", "phreeqc", "rds") - -init_cell <- list(H = 1.476571028625e-01, - O = 7.392297218936e-02, - Charge = -1.765225732724e-18, - `C(-4)` = 2.477908970828e-21, - `C(4)` = 2.647623016916e-06, - Ca = 2.889623169138e-05, - Cl = 4.292806181039e-04, - `Fe(2)` =1.908142472666e-07, - `Fe(3)` =3.173306589931e-12, - `H(0)` =2.675642675119e-15, - K = 2.530134809667e-06, - Mg =2.313806319294e-05, - Na =3.674633059628e-04, - `S(-2)` = 8.589766637180e-15, - `S(2)` = 1.205284362720e-19, - `S(4)` = 9.108958772790e-18, - `S(6)` = 2.198092329098e-05, - Sr = 6.012080128154e-07, - `U(4)` = 1.039668623852e-14, - `U(5)` = 1.208394829796e-15, - `U(6)` = 2.976409147150e-12) - -grid <- list( - n_cells = c(n, m), - s_cells = c(1, 1), - type = "scratch" +diffusion_setup <- list( + boundaries = list( + "N" = bound_def + ), + alpha_x = 1e-6, + alpha_y = 1e-6 ) -################################################################## -## Section 2 ## -## Diffusion parameters and boundary conditions ## -################################################################## - -vecinj_diffu <- list( - list(H = 0.147659686316291, - O = 0.0739242798146046, - Charge = 7.46361643222701e-20, - `C(-4)` = 2.92438561098248e-21, - `C(4)` = 2.65160558871092e-06, - Ca = 2.89001071336443e-05, - Cl = 0.000429291158114428, - `Fe(2)` = 1.90823391198114e-07, - `Fe(3)` = 3.10832423034763e-12, - `H(0)` = 2.7888235127385e-15, - K = 2.5301787e-06, - Mg = 2.31391999937907e-05, - Na = 0.00036746969, - `S(-2)` = 1.01376078438546e-14, - `S(2)` = 1.42247026981542e-19, - `S(4)` = 9.49422092568557e-18, - `S(6)` = 2.19812504654191e-05, - Sr = 6.01218519999999e-07, - `U(4)` = 4.82255946569383e-12, - `U(5)` = 5.49050615347901e-13, - `U(6)` = 1.32462838991902e-09) -) - -vecinj <- do.call(rbind.data.frame, vecinj_diffu) -names(vecinj) <- grid$props - -## diffusion coefficients -alpha_diffu <- c(H = 1E-6, O = 1E-6, Charge = 1E-6, `C(-4)` = 1E-6, - `C(4)` = 1E-6, Ca = 1E-6, Cl = 1E-6, `Fe(2)` = 1E-6, - `Fe(3)` = 1E-6, `H(0)` = 1E-6, K = 1E-6, Mg = 1E-6, - Na = 1E-6, `S(-2)` = 1E-6, `S(2)` = 1E-6, - `S(4)` = 1E-6, `S(6)` = 1E-6, Sr = 1E-6, - `U(4)` = 1E-6, `U(5)` = 1E-6, `U(6)` = 1E-6) - -## list of boundary conditions/inner nodes - -## vecinj_inner <- list( -## list(1,1,1) -## ) - -boundary <- list( - "N" = rep(1, n), - "E" = rep(0, n), - "S" = rep(0, n), - "W" = rep(0, n) -) - -diffu_list <- names(alpha_diffu) - -vecinj <- do.call(rbind.data.frame, vecinj_diffu) -names(vecinj) <- names(init_cell) - -diffusion <- list( - init = as.data.frame(init_cell, check.names = FALSE), - vecinj = vecinj, -# vecinj_inner = vecinj_inner, - vecinj_index = boundary, - alpha = alpha_diffu -) - -################################################################# -## Section 3 ## -## Chemistry module (Phreeqc) ## -################################################################# - - -chemistry <- list( - database = database, - input_script = input_script -) - -################################################################# -## Section 4 ## -## Putting all those things together ## -################################################################# - - -iterations <- 10 -dt <- 200 +chemistry_setup <- list() +# Define a setup list for simulation configuration setup <- list( - grid = grid, - diffusion = diffusion, - chemistry = chemistry, - iterations = iterations, - timesteps = rep(dt, iterations), - store_result = TRUE -) + Grid = grid_setup, # Parameters related to the grid structure + Diffusion = diffusion_setup, # Parameters related to the diffusion process + Chemistry = chemistry_setup # Parameters related to the chemistry process +) \ No newline at end of file diff --git a/bench/surfex/ex_rt.R b/bench/surfex/ex_rt.R new file mode 100644 index 000000000..dff0b5c1a --- /dev/null +++ b/bench/surfex/ex_rt.R @@ -0,0 +1,7 @@ +iterations <- 10 +dt <- 200 + +list( + timesteps = rep(dt, iterations), + store_result = TRUE +) \ No newline at end of file diff --git a/bench/surfex/surfex.R b/bench/surfex/surfex.R index 27409b97f..df05aab70 100644 --- a/bench/surfex/surfex.R +++ b/bench/surfex/surfex.R @@ -1,141 +1,37 @@ -## Time-stamp: "Last modified 2023-08-02 13:59:44 mluebke" +rows <- 1000 +cols <- 1000 -database <- normalizePath("../share/poet/bench/surfex/SMILE_2021_11_01_TH.dat") -input_script <- normalizePath("../share/poet/bench/surfex/SurfExBase.pqi") +grid_def <- matrix(1, nrow = rows, ncol = cols) -cat(paste(":: R This is a test 1\n")) +# Define grid configuration for POET model +grid_setup <- list( + pqc_in_file = "./SurfExBase.pqi", + pqc_db_file = "./SMILE_2021_11_01_TH.dat", # Path to the database file for Phreeqc + grid_def = grid_def, # Definition of the grid, containing IDs according to the Phreeqc input script + grid_size = c(rows, cols) / 10, # Size of the grid in meters + constant_cells = c() # IDs of cells with constant concentration +) -################################################################# -## Section 1 ## -## Grid initialization ## -################################################################# +bound_def <- list( + "type" = rep("constant", cols), + "sol_id" = rep(2, cols), + "cell" = seq(1, cols) +) -n <- 1000 -m <- 1000 - -types <- c("scratch", "phreeqc", "rds") - -init_cell <- list(H = 1.476571028625e-01, - O = 7.392297218936e-02, - Charge = -1.765225732724e-18, - `C(-4)` = 2.477908970828e-21, - `C(4)` = 2.647623016916e-06, - Ca = 2.889623169138e-05, - Cl = 4.292806181039e-04, - `Fe(2)` =1.908142472666e-07, - `Fe(3)` =3.173306589931e-12, - `H(0)` =2.675642675119e-15, - K = 2.530134809667e-06, - Mg =2.313806319294e-05, - Na =3.674633059628e-04, - `S(-2)` = 8.589766637180e-15, - `S(2)` = 1.205284362720e-19, - `S(4)` = 9.108958772790e-18, - `S(6)` = 2.198092329098e-05, - Sr = 6.012080128154e-07, - `U(4)` = 1.039668623852e-14, - `U(5)` = 1.208394829796e-15, - `U(6)` = 2.976409147150e-12) - -grid <- list( - n_cells = c(n, m), - s_cells = c(n/10, m/10), - type = "scratch" +diffusion_setup <- list( + boundaries = list( + "N" = bound_def + ), + alpha_x = 1e-6, + alpha_y = 1e-6 ) -################################################################## -## Section 2 ## -## Diffusion parameters and boundary conditions ## -################################################################## - -vecinj_diffu <- list( - list(H = 0.147659686316291, - O = 0.0739242798146046, - Charge = 7.46361643222701e-20, - `C(-4)` = 2.92438561098248e-21, - `C(4)` = 2.65160558871092e-06, - Ca = 2.89001071336443e-05, - Cl = 0.000429291158114428, - `Fe(2)` = 1.90823391198114e-07, - `Fe(3)` = 3.10832423034763e-12, - `H(0)` = 2.7888235127385e-15, - K = 2.5301787e-06, - Mg = 2.31391999937907e-05, - Na = 0.00036746969, - `S(-2)` = 1.01376078438546e-14, - `S(2)` = 1.42247026981542e-19, - `S(4)` = 9.49422092568557e-18, - `S(6)` = 2.19812504654191e-05, - Sr = 6.01218519999999e-07, - `U(4)` = 4.82255946569383e-12, - `U(5)` = 5.49050615347901e-13, - `U(6)` = 1.32462838991902e-09) -) - -vecinj <- do.call(rbind.data.frame, vecinj_diffu) -names(vecinj) <- grid$props - -## diffusion coefficients -alpha_diffu <- c(H = 1E-6, O = 1E-6, Charge = 1E-6, `C(-4)` = 1E-6, - `C(4)` = 1E-6, Ca = 1E-6, Cl = 1E-6, `Fe(2)` = 1E-6, - `Fe(3)` = 1E-6, `H(0)` = 1E-6, K = 1E-6, Mg = 1E-6, - Na = 1E-6, `S(-2)` = 1E-6, `S(2)` = 1E-6, - `S(4)` = 1E-6, `S(6)` = 1E-6, Sr = 1E-6, - `U(4)` = 1E-6, `U(5)` = 1E-6, `U(6)` = 1E-6) - -## list of boundary conditions/inner nodes - -## vecinj_inner <- list( -## list(1,1,1) -## ) - -boundary <- list( - "N" = rep(1, n), - "E" = rep(0, n), - "S" = rep(0, n), - "W" = rep(0, n) -) - -diffu_list <- names(alpha_diffu) - -vecinj <- do.call(rbind.data.frame, vecinj_diffu) -names(vecinj) <- names(init_cell) - -diffusion <- list( - init = as.data.frame(init_cell, check.names = FALSE), - vecinj = vecinj, -# vecinj_inner = vecinj_inner, - vecinj_index = boundary, - alpha = alpha_diffu -) - -################################################################# -## Section 3 ## -## Chemistry module (Phreeqc) ## -################################################################# - - -chemistry <- list( - database = database, - input_script = input_script -) - -################################################################# -## Section 4 ## -## Putting all those things together ## -################################################################# - - -iterations <- 100 -dt <- 200 +chemistry_setup <- list() +# Define a setup list for simulation configuration setup <- list( - grid = grid, - diffusion = diffusion, - chemistry = chemistry, - iterations = iterations, - timesteps = rep(dt, iterations), - store_result = TRUE, - out_save = c(5, iterations) -) + Grid = grid_setup, # Parameters related to the grid structure + Diffusion = diffusion_setup, # Parameters related to the diffusion process + Chemistry = chemistry_setup # Parameters related to the chemistry process +) \ No newline at end of file diff --git a/bench/surfex/surfex_rt.R b/bench/surfex/surfex_rt.R new file mode 100644 index 000000000..e429fd771 --- /dev/null +++ b/bench/surfex/surfex_rt.R @@ -0,0 +1,10 @@ +iterations <- 100 +dt <- 200 + +out_save <- seq(5, iterations, by = 5) + +list( + timesteps = rep(dt, iterations), + store_result = TRUE, + out_save = out_save +) \ No newline at end of file diff --git a/docs/CMakeLists.txt b/docs/CMakeLists.txt index 4781c2b87..fa9038767 100644 --- a/docs/CMakeLists.txt +++ b/docs/CMakeLists.txt @@ -14,7 +14,6 @@ if(DOXYGEN_FOUND) doxygen_add_docs(doxygen ${PROJECT_SOURCE_DIR}/src ${PROJECT_SOURCE_DIR}/README.md - ${PROJECT_SOURCE_DIR}/docs/Input_Scripts.md ${PROJECT_SOURCE_DIR}/docs/Output.md COMMENT "Generate html pages") endif() diff --git a/docs/Input_Scripts.md b/docs/Input_Scripts.md deleted file mode 100644 index 1c03de749..000000000 --- a/docs/Input_Scripts.md +++ /dev/null @@ -1,86 +0,0 @@ -# Input Scripts - -In the following the expected schemes of the input scripts is described. -Therefore, each section of the input script gets its own chapter. All sections -should return a `list` as results, which are concatenated to one setup list at -the end of the file. All values must have the same name in order to get parsed -by POET. - -## Grid initialization - -| name | type | description | -|----------------|----------------|-----------------------------------------------------------------------| -| `n_cells` | Numeric Vector | Number of cells in each direction | -| `s_cells` | Numeric Vector | Spatial resolution of grid in each direction | -| `type` | String | Type of initialization, can be set to *scratch*, *phreeqc* or *rds* | - -## Diffusion parameters - -| name | type | description | -|----------------|----------------------|-------------------------------------------| -| `init` | Named Numeric Vector | Initial state for each diffused species | -| `vecinj` | Data Frame | Defining all boundary conditions row wise | -| `vecinj_inner` | List of Triples | Inner boundaries | -| `vecinj_index` | List of 4 elements | Ghost nodes boundary conditions | -| `alpha` | Named Numeric Vector | Constant alpha for each species | - -### Remark on boundary conditions - -Each boundary condition should be defined in `vecinj` as a data frame, where one -row holds one boundary condition. - -To define inner (constant) boundary conditions, use a list of triples in -`vecinj_inner`, where each triples is defined by $(i,x,y)$. $i$ is defining the -boundary condition, referencing to the row in `vecinj`. $x$ and $y$ coordinates -then defining the position inside the grid. - -Ghost nodes are set by `vecinj_index` which is a list containing boundaries for -each celestial direction (**important**: named by `N, E, S, W`). Each direction -is a numeric vector, also representing a row index of the `vecinj` data frame -for each ghost node, starting at the left-most and upper cell respectively. By -setting the boundary condition to $0$, the ghost node is set as closed boundary. - -#### Example - -Suppose you have a `vecinj` data frame defining 2 boundary conditions and a grid -consisting of $10 \times 10$ grid cells. Grid cell $(1,1)$ should be set to the -first boundary condition and $(5,6)$ to the second. Also, all boundary -conditions for the ghost nodes should be closed. Except the southern boundary, -which should be set to the first boundary condition injection. The following -setup describes how to setup your initial script, where `n` and `m` are the -grids cell count for each direction ($n = m = 10$): - -```R -vecinj_inner <- list ( - l1 = c(1, 1, 1), - l2 = c(2, 5, 6) -) - -vecinj_index <- list( - "N" = rep(0, n), - "E" = rep(0, m), - "S" = rep(1, n), - "W" = rep(0, m) -) -``` - -## Chemistry parameters - -| name | type | description | -|----------------|--------------|----------------------------------------------------------------------------------| -| `database` | String | Path to the Phreeqc database | -| `input_script` | String | Path the the Phreeqc input script | -| `dht_species` | Named Vector | Indicates significant digits to use for each species for DHT rounding. | -| `pht_species` | Named Vector | Indicates significant digits to use for each species for Interpolation rounding. | - -## Final setup - -| name | type | description | -|----------------|----------------|------------------------------------------------------------| -| `grid` | List | Grid parameter list | -| `diffusion` | List | Diffusion parameter list | -| `chemistry` | List | Chemistry parameter list | -| `iterations` | Numeric Value | Count of iterations | -| `timesteps` | Numeric Vector | $\Delta t$ to use for specific iteration | -| `store_result` | Boolean | Indicates if results should be stored | -| `out_save` | Numeric Vector | *optional:* At which iteration the states should be stored | diff --git a/docs/Output.md b/docs/Output.md index 4044c84a6..82644d347 100644 --- a/docs/Output.md +++ b/docs/Output.md @@ -35,34 +35,50 @@ corresponding values can be found in `/timings.rds` and possible to read out within a R runtime with `readRDS("timings.rds")`. There you will find the following values: -| Value | Description | -|--------------------|----------------------------------------------------------------------------| -| simtime | time spent in whole simulation loop without any initialization and cleanup | -| simtime\_transport | measured time in *transport* subroutine | -| simtime\_chemistry | measured time in *chemistry* subroutine (actual parallelized part) | +| Value | Description | +| --------- | -------------------------------------------------------------------------- | +| simtime | time spent in whole simulation loop without any initialization and cleanup | +| chemistry | measured time in *chemistry* subroutine | +| diffusion | measured time in *diffusion* subroutine | -### chemistry subsetting +### Chemistry subsetting -If running parallel there are also measured timings which are subsets of -*simtime\_​chemistry*. +| Value | Description | +| ------------- | --------------------------------------------------------- | +| simtime | overall runtime of chemistry | +| loop | time spent in send/recv loop of master | +| sequential | sequential part of the master (e.g. shuffling field) | +| idle\_master | idling time of the master waiting for workers | +| idle\_worker | idling time (waiting for work from master) of the workers | +| phreeqc\_time | accumulated times for Phreeqc calls of every worker | -| Value | Description | -|-----------------------|-----------------------------------------------------------| -| chemistry\_loop | time spent in send/​recv loop of master | -| chemistry\_sequential | sequential part of master chemistry | -| idle\_master | idling time (waiting for any free worker) of the master | -| idle\_worker | idling time (waiting for work from master) of the workers | -| phreeqc\_time | accumulated times for Phreeqc calls of every worker | - -### DHT usage {#DHT-usage} +#### DHT usage If running in parallel and with activated DHT, two more timings and also some profiling about the DHT usage are given: | Value | Description | -|-----------------|---------------------------------------------------------| -| dht\_fill\_time | time to write data to DHT | -| dht\_get\_time | time to retreive data from DHT | -| dh\_hits | count of data points retrieved from DHT | -| dht\_miss | count of misses/count of data points written to DHT | +| --------------- | ------------------------------------------------------- | +| dht\_hits | count of data points retrieved from DHT | | dht\_evictions | count of data points evicted by another write operation | +| dht\_get\_time | time to retreive data from DHT | +| dht\_fill\_time | time to write data to DHT | + +#### Interpolation + +If using interpolation, the following values are given: + +| Value | Description | +| -------------- | --------------------------------------------------------------------- | +| interp\_w | time spent to write to PHT | +| interp\_r | time spent to read from DHT/PHT/Cache | +| interp\_g | time spent to gather results from DHT | +| interp\_fc | accumulated time spent in interpolation function call | +| interp\_calls | count of interpolations | +| interp\_cached | count of interpolation data sets, which where cached in the local map | + +### Diffusion subsetting + +| Value | Description | +| --------- | ------------------------------------------ | +| simtime | overall runtime of diffusion | diff --git a/ext/doctest b/ext/doctest deleted file mode 160000 index ae7a13539..000000000 --- a/ext/doctest +++ /dev/null @@ -1 +0,0 @@ -Subproject commit ae7a13539fb71f270b87eb2e874fbac80bc8dda2 diff --git a/ext/iphreeqc b/ext/iphreeqc new file mode 160000 index 000000000..749fdbc2e --- /dev/null +++ b/ext/iphreeqc @@ -0,0 +1 @@ +Subproject commit 749fdbc2e9478046bf3f270c70e5800637246712 diff --git a/ext/phreeqcrm b/ext/phreeqcrm deleted file mode 160000 index 6ed14c353..000000000 --- a/ext/phreeqcrm +++ /dev/null @@ -1 +0,0 @@ -Subproject commit 6ed14c35322a245e3a9776ef262c0ac0eba3b301 diff --git a/ext/tug b/ext/tug index 25855da6b..449647010 160000 --- a/ext/tug +++ b/ext/tug @@ -1 +1 @@ -Subproject commit 25855da6b2930559b542bbadb16299932332d6a3 +Subproject commit 449647010ab9cdf9e405139f360424a2b21ab3ab diff --git a/src/Base/Grid.cpp b/src/Base/Grid.cpp deleted file mode 100644 index e51b114ba..000000000 --- a/src/Base/Grid.cpp +++ /dev/null @@ -1,274 +0,0 @@ -/* -** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of -** Potsdam) -** -** Copyright (C) 2018-2023 Marco De Lucia, Max Luebke (GFZ Potsdam) -** -** POET is free software; you can redistribute it and/or modify it under the -** terms of the GNU General Public License as published by the Free Software -** Foundation; either version 2 of the License, or (at your option) any later -** version. -** -** POET is distributed in the hope that it will be useful, but WITHOUT ANY -** WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR -** A PARTICULAR PURPOSE. See the GNU General Public License for more details. -** -** You should have received a copy of the GNU General Public License along with -** this program; if not, write to the Free Software Foundation, Inc., 51 -** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. -*/ - -#include "Grid.hpp" -#include "SimParams.hpp" - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -using namespace poet; -using namespace Rcpp; - -enum { INIT_SCRATCH, INIT_PHREEQC, INIT_RDS }; - -static inline int8_t get_type_id(std::string type) { - if (type == "scratch") { - return INIT_SCRATCH; - } - if (type == "phreeqc") { - return INIT_PHREEQC; - } - if (type == "rds") { - return INIT_RDS; - } - return -1; -} - -void poet::Grid::InitModuleFromParams(const poet::GridParams &grid_args) { - assert(grid_args.n_cells.size() == grid_args.s_cells.size()); - - this->SetGridDimension(grid_args.n_cells.size()); - - std::copy(grid_args.n_cells.begin(), grid_args.n_cells.end(), - this->n_cells.begin()); - std::copy(grid_args.s_cells.begin(), grid_args.s_cells.end(), - this->grid_size.begin()); - - auto init_type = get_type_id(grid_args.type); - assert(init_type >= 0); - - switch (init_type) { - case INIT_SCRATCH: - this->SetPropNames(grid_args.props); - this->InitGridFromScratch(grid_args.init_df); - break; - case INIT_RDS: - this->SetPropNames(grid_args.props); - this->InitGridFromRDS(); - break; - case INIT_PHREEQC: - this->InitGridFromPhreeqc(); - } -} - -void poet::Grid::SetGridDimension(uint8_t dim) { - assert(dim < MAX_DIM + 1); - - this->dim = dim; -} - -void poet::Grid::SetGridCellCount(uint32_t n_x, uint32_t n_y, uint32_t n_z) { - assert(this->dim < MAX_DIM + 1); - - this->n_cells[0] = n_x; - this->n_cells[1] = n_y; - // this->n_cells[2] = n_z; -} - -void poet::Grid::SetGridSize(double s_x, double s_y, double s_z) { - assert(this->dim < MAX_DIM + 1); - - this->grid_size[0] = s_x; - this->grid_size[1] = s_y; - // this->grid_size[2] = s_z; -} - -void poet::Grid::SetPropNames(const std::vector &prop_names) { - this->prop_names = prop_names; -} - -void poet::Grid::PushbackModuleFlow(const std::string &input, - const std::string &output) { - FlowInputOutputInfo element = {input, output}; - this->flow_vec.push_back(element); -} - -void poet::Grid::PreModuleFieldCopy(uint32_t tick) { - FlowInputOutputInfo curr_element = this->flow_vec.at(tick); - - const std::string input_module_name = curr_element.input_field; - StateMemory *out_state = this->GetStatePointer(curr_element.output_field); - - std::vector &mod_props = out_state->props; - std::vector &mod_field = out_state->mem; - - // copy output of another module as input for another module - for (uint32_t i = 0; i < mod_props.size(); i++) { - try { - std::vector t_prop_vec = - this->GetSpeciesByName(mod_props[i], input_module_name); - - std::copy(t_prop_vec.begin(), t_prop_vec.end(), - mod_field.begin() + (i * this->GetTotalCellCount())); - } catch (...) { - // TODO: there might be something wrong ... - continue; - } - } -} - -Grid::Grid() { - this->n_cells.fill(0); - this->grid_size.fill(0); -}; - -Grid::~Grid() { - for (auto &[key, val] : this->state_register) { - delete val; - } -} - -void poet::Grid::InitGridFromScratch(const Rcpp::DataFrame &init_cell) { - const uint32_t vec_size = this->GetTotalCellCount(); - - StateMemory *curr_state = - this->RegisterState(poet::GRID_MODULE_NAME, this->prop_names); - curr_state->props = this->prop_names; - - std::vector &curr_field = curr_state->mem; - - for (auto const &prop : this->prop_names) { - // size of the vector shall be 1 - double prop_val = ((Rcpp::NumericVector)init_cell[prop.c_str()])[0]; - std::vector prop_vec(vec_size, prop_val); - - curr_field.insert(curr_field.end(), prop_vec.begin(), prop_vec.end()); - } - - this->grid_init = curr_state; -} -void poet::Grid::InitGridFromRDS() { - // TODO -} -void poet::Grid::InitGridFromPhreeqc() { - // TODO -} - -poet::StateMemory *Grid::RegisterState(std::string module_name, - std::vector props) { - poet::StateMemory *mem = new poet::StateMemory; - - mem->props = props; - - const auto [it, success] = this->state_register.insert({module_name, mem}); - - if (!success) { - delete mem; - throw std::bad_alloc(); - } - - return mem; -} - -poet::StateMemory *Grid::GetStatePointer(std::string module_name) { - const auto it = this->state_register.find(module_name); - - if (it == this->state_register.end()) { - throw std::range_error( - std::string("Requested module " + module_name + " was not found")); - } - - return it->second; -} - -auto Grid::GetGridSize(uint8_t direction) const -> uint32_t { - return this->grid_size.at(direction); -} - -auto Grid::GetGridDimension() const -> uint8_t { return this->dim; } - -auto Grid::GetTotalCellCount() const -> uint32_t { - uint32_t sum = 1; - for (auto const &dim : this->n_cells) { - sum *= (dim != 0 ? dim : 1); - } - - return sum; -} - -auto Grid::GetGridCellsCount(uint8_t direction) const -> uint32_t { - return this->n_cells.at(direction); -} - -auto poet::Grid::GetInitialGrid() const -> StateMemory * { - return this->grid_init; -} - -auto Grid::GetSpeciesCount() const -> uint32_t { - return this->prop_names.size(); -} - -auto Grid::GetPropNames() const -> std::vector { - return this->prop_names; -} - -auto poet::Grid::GetSpeciesByName(std::string name, - std::string module_name) const - -> std::vector { - - auto const it = this->state_register.find(module_name); - - if (it == this->state_register.end()) { - throw std::range_error( - std::string("Requested module " + module_name + " was not found")); - } - - poet::StateMemory const *module_memory = it->second; - - auto const &prop_vector = module_memory->props; - auto const find_res = std::find(prop_vector.begin(), prop_vector.end(), name); - if (find_res == prop_vector.end()) { - throw std::range_error(std::string( - "Species " + name + " was not found for module " + module_name)); - } - uint32_t prop_index = find_res - prop_vector.begin(); - - uint32_t begin_vec = prop_index * this->GetTotalCellCount(), - end_vec = ((prop_index + 1) * this->GetTotalCellCount()); - - return std::vector(module_memory->mem.begin() + begin_vec, - module_memory->mem.begin() + end_vec); -} - -// HACK: Helper function -void poet::Grid::WriteFieldsToR(RInside &R) { - - for (auto const &elem : this->state_register) { - std::string field_name = elem.first; - StateMemory *field = elem.second; - - const uint32_t n_cells_per_prop = field->mem.size() / field->props.size(); - - R["TMP"] = Rcpp::wrap(field->mem); - R["TMP_PROPS"] = Rcpp::wrap(field->props); - R.parseEval(std::string( - "mysetup$" + field_name + " <- setNames(data.frame(matrix(TMP, nrow=" + - std::to_string(n_cells_per_prop) + ")), TMP_PROPS)")); - } -} diff --git a/src/Base/Grid.hpp b/src/Base/Grid.hpp deleted file mode 100644 index deec5c04a..000000000 --- a/src/Base/Grid.hpp +++ /dev/null @@ -1,117 +0,0 @@ -/* -** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of -** Potsdam) -** -** Copyright (C) 2018-2023 Marco De Lucia, Max Luebke (GFZ Potsdam) -** -** POET is free software; you can redistribute it and/or modify it under the -** terms of the GNU General Public License as published by the Free Software -** Foundation; either version 2 of the License, or (at your option) any later -** version. -** -** POET is distributed in the hope that it will be useful, but WITHOUT ANY -** WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR -** A PARTICULAR PURPOSE. See the GNU General Public License for more details. -** -** You should have received a copy of the GNU General Public License along with -** this program; if not, write to the Free Software Foundation, Inc., 51 -** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. -*/ - -#ifndef GRID_H -#define GRID_H - -#include "SimParams.hpp" - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#define MAX_DIM 2 - -namespace poet { - -enum { GRID_X_DIR, GRID_Y_DIR, GRID_Z_DIR }; - -using StateMemory = struct s_StateMemory { - std::vector mem; - std::vector props; -}; - -using FlowInputOutputInfo = struct inOut_info { - std::string input_field; - std::string output_field; -}; - -constexpr const char *GRID_MODULE_NAME = "grid_init"; - -/** - * @brief Class describing the grid - * - * Providing methods to shuffle and unshuffle grid (for the master) as also to - * import and export a work package (for worker). - * - * @todo find better abstraction - * - */ -class Grid { - -public: - Grid(); - - ~Grid(); - - void InitModuleFromParams(const poet::GridParams &grid_args); - - void SetGridDimension(uint8_t dim); - void SetGridCellCount(uint32_t n_x, uint32_t n_y = 0, uint32_t n_z = 0); - void SetGridSize(double s_x, double s_y = 0., double s_z = 0.); - void SetPropNames(const std::vector &prop_names); - - void PushbackModuleFlow(const std::string &input, const std::string &output); - void PreModuleFieldCopy(uint32_t tick); - - void InitGridFromScratch(const Rcpp::DataFrame &init_cell); - void InitGridFromRDS(); - void InitGridFromPhreeqc(); - - auto GetGridDimension() const -> uint8_t; - auto GetTotalCellCount() const -> uint32_t; - auto GetGridCellsCount(uint8_t direction) const -> uint32_t; - auto GetGridSize(uint8_t direction) const -> uint32_t; - - auto RegisterState(std::string module_name, std::vector props) - -> StateMemory *; - auto GetStatePointer(std::string module_name) -> StateMemory *; - - auto GetInitialGrid() const -> StateMemory *; - - auto GetSpeciesCount() const -> uint32_t; - auto GetPropNames() const -> std::vector; - - auto GetSpeciesByName(std::string name, - std::string module_name = poet::GRID_MODULE_NAME) const - -> std::vector; - - void WriteFieldsToR(RInside &R); - -private: - std::vector flow_vec; - - std::uint8_t dim = 0; - std::array grid_size; - std::array n_cells; - - std::map state_register; - StateMemory *grid_init = std::nullptr_t(); - - std::vector prop_names; -}; -} // namespace poet -#endif // GRID_H diff --git a/src/Base/RInsidePOET.hpp b/src/Base/RInsidePOET.hpp index 6467e5760..2897fc5a8 100644 --- a/src/Base/RInsidePOET.hpp +++ b/src/Base/RInsidePOET.hpp @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -42,6 +43,13 @@ public: } } + RHookFunction(SEXP f) { + try { + this->func = Rcpp::Function(f); + } catch (const std::exception &e) { + } + } + template T operator()(Args... args) const { if (func.has_value()) { return (Rcpp::as(this->func.value()(args...))); @@ -50,8 +58,17 @@ public: } } + RHookFunction &operator=(const RHookFunction &rhs) { + this->func = rhs.func; + return *this; + } + + RHookFunction(const RHookFunction &rhs) { this->func = rhs.func; } + bool isValid() const { return this->func.has_value(); } + SEXP asSEXP() const { return Rcpp::as(this->func.value()); } + private: std::optional func; }; diff --git a/src/Base/SimParams.cpp b/src/Base/SimParams.cpp deleted file mode 100644 index f711e796a..000000000 --- a/src/Base/SimParams.cpp +++ /dev/null @@ -1,266 +0,0 @@ -// Time-stamp: "Last modified 2023-08-15 12:12:36 delucia" - -/* -** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of -** Potsdam) -** -** Copyright (C) 2018-2023 Marco De Lucia, Max Luebke (GFZ Potsdam) -** -** POET is free software; you can redistribute it and/or modify it under the -** terms of the GNU General Public License as published by the Free Software -** Foundation; either version 2 of the License, or (at your option) any later -** version. -** -** POET is distributed in the hope that it will be useful, but WITHOUT ANY -** WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR -** A PARTICULAR PURPOSE. See the GNU General Public License for more details. -** -** You should have received a copy of the GNU General Public License along with -** this program; if not, write to the Free Software Foundation, Inc., 51 -** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. -*/ - -#include "SimParams.hpp" - -#include "../Chemistry/enums.hpp" - -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -using namespace poet; -using namespace std; -using namespace Rcpp; - -poet::GridParams::s_GridParams(RInside &R) { - auto tmp_n_cells = - Rcpp::as>(R.parseEval("mysetup$grid$n_cells")); - assert(tmp_n_cells.size() < 3); - - this->dim = tmp_n_cells.size(); - - std::copy(tmp_n_cells.begin(), tmp_n_cells.end(), this->n_cells.begin()); - - auto tmp_s_cells = - Rcpp::as>(R.parseEval("mysetup$grid$s_cells")); - - assert(tmp_s_cells.size() == this->dim); - - std::copy(tmp_s_cells.begin(), tmp_s_cells.end(), this->s_cells.begin()); - - this->total_n = - (dim == 1 ? this->n_cells[0] : this->n_cells[0] * this->n_cells[1]); - - this->type = Rcpp::as(R.parseEval("mysetup$grid$type")); -} - -poet::DiffusionParams::s_DiffusionParams(RInside &R) { - this->initial_t = - Rcpp::as(R.parseEval("mysetup$diffusion$init")); - this->alpha = - Rcpp::as(R.parseEval("mysetup$diffusion$alpha")); - if (Rcpp::as( - R.parseEval("'vecinj_inner' %in% names(mysetup$diffusion)"))) { - this->vecinj_inner = - Rcpp::as(R.parseEval("mysetup$diffusion$vecinj_inner")); - } - this->vecinj = - Rcpp::as(R.parseEval("mysetup$diffusion$vecinj")); - this->vecinj_index = - Rcpp::as(R.parseEval("mysetup$diffusion$vecinj_index")); -} - -void poet::ChemistryParams::initFromR(RInsidePOET &R) { - this->database_path = - Rcpp::as(R.parseEval("mysetup$chemistry$database")); - this->input_script = - Rcpp::as(R.parseEval("mysetup$chemistry$input_script")); - - if (R.checkIfExists("dht_species", "mysetup$chemistry")) { - this->dht_signifs = Rcpp::as>( - R.parseEval(("mysetup$chemistry$dht_species"))); - } - - if (R.checkIfExists("pht_species", "mysetup$chemistry")) { - this->pht_signifs = Rcpp::as>( - R.parseEval(("mysetup$chemistry$pht_species"))); - } - this->hooks.dht_fill = - RHookFunction(R, "mysetup$chemistry$hooks$dht_fill"); - this->hooks.dht_fuzz = - RHookFunction>(R, "mysetup$chemistry$hooks$dht_fuzz"); - this->hooks.interp_pre = RHookFunction>( - R, "mysetup$chemistry$hooks$interp_pre_func"); - this->hooks.interp_post = - RHookFunction(R, "mysetup$chemistry$hooks$interp_post_func"); -} - -SimParams::SimParams(int world_rank_, int world_size_) { - this->simparams.world_rank = world_rank_; - this->simparams.world_size = world_size_; -} - -int SimParams::parseFromCmdl(char *argv[], RInsidePOET &R) { - // initialize argh object - argh::parser cmdl(argv); - - // if user asked for help - if (cmdl[{"help", "h"}]) { - if (simparams.world_rank == 0) { - MSG("Todo"); - MSG("See README.md for further information."); - } - return poet::PARSER_HELP; - } - // if positional arguments are missing - else if (!cmdl(2)) { - if (simparams.world_rank == 0) { - ERRMSG("Kin needs 2 positional arguments: "); - ERRMSG("1) the R script defining your simulation and"); - ERRMSG("2) the directory prefix where to save results and profiling"); - } - return poet::PARSER_ERROR; - } - - // collect all parameters which are not known, print them to stderr and return - // with PARSER_ERROR - std::list optionsError = validateOptions(cmdl); - if (!optionsError.empty()) { - if (simparams.world_rank == 0) { - ERRMSG("Unrecognized option(s):\n"); - for (auto option : optionsError) { - ERRMSG(std::string(option)); - } - ERRMSG("Make sure to use available options. Exiting!"); - } - return poet::PARSER_ERROR; - } - - simparams.print_progressbar = cmdl[{"P", "progress"}]; - - // simparams.print_progressbar = cmdl[{"P", "progress"}]; - - /* Parse DHT arguments */ - chem_params.use_dht = cmdl["dht"]; - chem_params.use_interp = cmdl["interp"]; - // cout << "CPP: DHT is " << ( dht_enabled ? "ON" : "OFF" ) << '\n'; - - cmdl("dht-size", DHT_SIZE_PER_PROCESS_MB) >> chem_params.dht_size; - // cout << "CPP: DHT size per process (Byte) = " << dht_size_per_process << - // endl; - - cmdl("dht-snaps", 0) >> chem_params.dht_snaps; - - cmdl("dht-file") >> chem_params.dht_file; - /*Parse work package size*/ - cmdl("work-package-size", WORK_PACKAGE_SIZE_DEFAULT) >> simparams.wp_size; - - cmdl("interp-size", 100) >> chem_params.pht_size; - cmdl("interp-min", 5) >> chem_params.interp_min_entries; - cmdl("interp-bucket-entries", 20) >> chem_params.pht_max_entries; - - /*Parse output options*/ - simparams.store_result = !cmdl["ignore-result"]; - - /*Parse work package size*/ - cmdl("work-package-size", WORK_PACKAGE_SIZE_DEFAULT) >> simparams.wp_size; - - chem_params.use_interp = cmdl["interp"]; - cmdl("interp-size", 100) >> chem_params.pht_size; - cmdl("interp-min", 5) >> chem_params.interp_min_entries; - cmdl("interp-bucket-entries", 20) >> chem_params.pht_max_entries; - - /*Parse output options*/ - simparams.store_result = !cmdl["ignore-result"]; - - if (simparams.world_rank == 0) { - MSG("Complete results storage is " + BOOL_PRINT(simparams.store_result)); - MSG("Work Package Size: " + std::to_string(simparams.wp_size)); - MSG("DHT is " + BOOL_PRINT(chem_params.use_dht)); - - if (chem_params.use_dht) { - MSG("DHT strategy is " + std::to_string(simparams.dht_strategy)); - // MDL: these should be outdated (?) - // MSG("DHT key default digits (ignored if 'signif_vector' is " - // "defined) = " - // << simparams.dht_significant_digits); - // MSG("DHT logarithm before rounding: " - // << (simparams.dht_log ? "ON" : "OFF")); - MSG("DHT size per process (Megabyte) = " + - std::to_string(chem_params.dht_size)); - MSG("DHT save snapshots is " + BOOL_PRINT(chem_params.dht_snaps)); - MSG("DHT load file is " + chem_params.dht_file); - } - - if (chem_params.use_interp) { - MSG("PHT interpolation enabled: " + BOOL_PRINT(chem_params.use_interp)); - MSG("PHT interp-size = " + std::to_string(chem_params.pht_size)); - MSG("PHT interp-min = " + - std::to_string(chem_params.interp_min_entries)); - MSG("PHT interp-bucket-entries = " + - std::to_string(chem_params.pht_max_entries)); - } - } - - cmdl(1) >> filesim; - cmdl(2) >> out_dir; - - chem_params.dht_outdir = out_dir; - - /* distribute information to R runtime */ - // if local_rank == 0 then master else worker - R["local_rank"] = simparams.world_rank; - // assign a char* (string) to 'filesim' - R["filesim"] = wrap(filesim); - // assign a char* (string) to 'fileout' - R["fileout"] = wrap(out_dir); - // pass the boolean "store_result" to the R process - R["store_result"] = simparams.store_result; - // worker count - R["n_procs"] = simparams.world_size - 1; - // work package size - R["work_package_size"] = simparams.wp_size; - // dht enabled? - R["dht_enabled"] = chem_params.use_dht; - // log before rounding? - R["dht_log"] = simparams.dht_log; - - // eval the init string, ignoring any returns - R.parseEvalQ("source(filesim)"); - R.parseEvalQ("mysetup <- setup"); - - this->chem_params.initFromR(R); - - return poet::PARSER_OK; -} - -void SimParams::initVectorParams(RInside &R) {} - -std::list SimParams::validateOptions(argh::parser cmdl) { - /* store all unknown parameters here */ - std::list retList; - - /* loop over all flags and compare to given flaglist*/ - for (auto &flag : cmdl.flags()) { - if (!(flaglist.find(flag) != flaglist.end())) - retList.push_back(flag); - } - - /* and loop also over params and compare to given paramlist */ - for (auto ¶m : cmdl.params()) { - if (!(paramlist.find(param.first) != paramlist.end())) - retList.push_back(param.first); - } - - return retList; -} diff --git a/src/Base/SimParams.hpp b/src/Base/SimParams.hpp deleted file mode 100644 index 523b5b105..000000000 --- a/src/Base/SimParams.hpp +++ /dev/null @@ -1,275 +0,0 @@ -/* -** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of -** Potsdam) -** -** Copyright (C) 2018-2022 Marco De Lucia, Max Luebke (GFZ Potsdam) -** -** POET is free software; you can redistribute it and/or modify it under the -** terms of the GNU General Public License as published by the Free Software -** Foundation; either version 2 of the License, or (at your option) any later -** version. -** -** POET is distributed in the hope that it will be useful, but WITHOUT ANY -** WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR -** A PARTICULAR PURPOSE. See the GNU General Public License for more details. -** -** You should have received a copy of the GNU General Public License along with -** this program; if not, write to the Free Software Foundation, Inc., 51 -** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. -*/ - -#ifndef PARSER_H -#define PARSER_H - -#include "../DataStructures/DataStructures.hpp" -#include "Macros.hpp" -#include "RInsidePOET.hpp" - -#include "argh.hpp" // Argument handler https://github.com/adishavit/argh - -#include -#include -#include -#include -#include -#include - -#include -#include -// BSD-licenced - -/** Standard DHT Size. Defaults to 1 GB (1000 MB) */ -constexpr uint32_t DHT_SIZE_PER_PROCESS_MB = 1.5E3; -/** Standard work package size */ -#define WORK_PACKAGE_SIZE_DEFAULT 32 - -namespace poet { - -enum { PARSER_OK, PARSER_ERROR, PARSER_HELP }; - -/** - * @brief Defining all simulation parameters - * - */ -typedef struct { - /** Count of processes in MPI_COMM_WORLD */ - int world_size; - /** rank of proces in MPI_COMM_WORLD */ - int world_rank; - /** indicates if DHT should be used */ - bool dht_enabled; - /** apply logarithm to key before rounding */ - bool dht_log; - /** indicates if timestep dt differs between iterations */ - bool dt_differ; - /** Indicates, when a DHT snapshot should be written */ - int dht_snaps; - /** not implemented: How a DHT is distributed over processes */ - int dht_strategy; - /** Size of DHt per process in byter */ - unsigned int dht_size_per_process; - /** Default significant digit for rounding */ - int dht_significant_digits; - /** Default work package size */ - unsigned int wp_size; - /** indicates if resulting grid should be stored after every iteration */ - bool store_result; - /** indicating whether the progress bar during chemistry simulation should be - * printed or not */ - bool print_progressbar; - - bool interp_enabled; -} t_simparams; - -using GridParams = struct s_GridParams { - std::array n_cells; - std::array s_cells; - std::uint8_t dim; - std::uint32_t total_n; - - std::string type; - - Rcpp::DataFrame init_df; - std::string input_script; - std::string database_path; - - std::vector props; - - s_GridParams(RInside &R); -}; - -using DiffusionParams = struct s_DiffusionParams { - Rcpp::DataFrame initial_t; - - Rcpp::NumericVector alpha; - Rcpp::List vecinj_inner; - - Rcpp::DataFrame vecinj; - Rcpp::DataFrame vecinj_index; - - s_DiffusionParams(RInside &R); -}; - -struct ChemistryParams { - std::string database_path; - std::string input_script; - - bool use_dht; - std::uint64_t dht_size; - int dht_snaps; - std::string dht_file; - std::string dht_outdir; - NamedVector dht_signifs; - - bool use_interp; - std::uint64_t pht_size; - std::uint32_t pht_max_entries; - std::uint32_t interp_min_entries; - NamedVector pht_signifs; - - struct Chem_Hook_Functions { - RHookFunction dht_fill; - RHookFunction> dht_fuzz; - RHookFunction> interp_pre; - RHookFunction interp_post; - } hooks; - - void initFromR(RInsidePOET &R); -}; - -/** - * @brief Reads information from program arguments and R runtime - * - * Providing functions to initialize parameters of the simulation using command - * line parameters and parameters from the R runtime. This class will also parse - * arguments from the commandline and decides if argument is known or unknown. - * - * Stores and distribute current simulation parameters at any time. - * - */ -class SimParams { -public: - /** - * @brief Construct a new SimParams object - * - * With all given parameters a new instance of this class will be created. - * - * @param world_rank Rank of process inside MPI_COMM_WORLD - * @param world_size Size of communicator MPI_COMM_WORLD - */ - SimParams(int world_rank, int world_size); - - /** - * @brief Parse program arguments - * - * This is done by the argh.h library. - * - * First, the function will check if there is a flag 'help' or 'h'. If this is - * the case a help message is printed and the function will return with - * PARSER_HELP. - * - * Second, if there are not 2 positional arguments an error will be printed to - * stderr and the function returns with PARSER_ERROR. - * - * Then all given program parameters and flags will be read and checked, if - * there are known by validateOptions. A list of all unknown options might be - * returned, printed out and the function will return with PARSER_ERROR. - * Oterhwise the function continuos. - * - * Now all program arguments will be stored inside t_simparams struct, printed - * out and the function returns with PARSER_OK. - * - * Also, all parsed agruments are distributed to the R runtime. - * - * @param argv Argument value of the program - * @param R Instantiated R runtime - * @return int Returns with 0 if no error occured, otherwise value less than 0 - * is returned. - */ - int parseFromCmdl(char *argv[], RInsidePOET &R); - - /** - * @brief Init std::vector values - * - * This will initialize dht_signif_vector and dht_prop_type_vector internally - * depending on whether vectors are defined by R-Simulation file or not. - * 'init_chemistry' must be run beforehand. - * - * @param R R runtime - */ - void initVectorParams(RInside &R); - - /** - * @brief Get the numerical params struct - * - * Returns a struct which contains all numerical or boolean simulation - * parameters. - * - * @return t_simparams Parameter struct - */ - auto getNumParams() const { return this->simparams; }; - - /** - * @brief Get the DHT_Signif_Vector - * - * Returns a vector indicating which significant values are used for each - * variable of a grid cell. - * - * @return std::vector Vector of integers containing information about - * significant digits - */ - auto getDHTSignifVector() const { return this->dht_signif_vector; }; - - auto getPHTSignifVector() const { return this->pht_signif_vector; }; - - auto getPHTBucketSize() const { return this->pht_bucket_size; }; - auto getPHTMinEntriesNeeded() const { return this->pht_min_entries_needed; }; - - /** - * @brief Get the filesim name - * - * Returns a string containing the absolute path to a R file defining the - * simulation. - * - * @return std::string Absolute path to R file - */ - auto getFilesim() const { return this->filesim; }; - - /** - * @brief Get the output directory - * - * Returns the name of an absolute path where all output files should be - * stored. - * - * @return std::string Absolute path to output path - */ - auto getOutDir() const { return this->out_dir; }; - - const auto &getChemParams() const { return chem_params; } - -private: - std::list validateOptions(argh::parser cmdl); - - const std::set flaglist{"ignore-result", "dht", "P", "progress", - "interp"}; - const std::set paramlist{ - "work-package-size", "dht-strategy", - "dht-size", "dht-snaps", - "dht-file", "interp-size", - "interp-min", "interp-bucket-entries"}; - - t_simparams simparams; - - std::vector dht_signif_vector; - std::vector pht_signif_vector; - - uint32_t pht_bucket_size; - uint32_t pht_min_entries_needed; - - std::string filesim; - std::string out_dir; - - ChemistryParams chem_params; -}; -} // namespace poet -#endif // PARSER_H diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 1494c0849..5d77d60e0 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,49 +1,85 @@ -add_library(poetlib - Base/Grid.cpp - Base/SimParams.cpp - Chemistry/ChemistryModule.cpp - Chemistry/MasterFunctions.cpp - Chemistry/WorkerFunctions.cpp - Chemistry/SurrogateModels/DHT_Wrapper.cpp - Chemistry/SurrogateModels/DHT.c - Chemistry/SurrogateModels/HashFunctions.cpp - Chemistry/SurrogateModels/InterpolationModule.cpp - Chemistry/SurrogateModels/ProximityHashTable.cpp +add_library(POETLib) + +target_sources(POETLib + PRIVATE + Init/InitialList.cpp + Init/GridInit.cpp + Init/DiffusionInit.cpp + Init/ChemistryInit.cpp DataStructures/Field.cpp Transport/DiffusionModule.cpp + Chemistry/ChemistryModule.cpp + Chemistry/MasterFunctions.cpp + Chemistry/WorkerFunctions.cpp + Chemistry/SurrogateModels/DHT_Wrapper.cpp + Chemistry/SurrogateModels/DHT.c + Chemistry/SurrogateModels/HashFunctions.cpp + Chemistry/SurrogateModels/InterpolationModule.cpp + Chemistry/SurrogateModels/ProximityHashTable.cpp ) -target_link_libraries(poetlib PUBLIC - MPI::MPI_C - ${MATH_LIBRARY} - RRuntime - PhreeqcRM - tug +target_include_directories(POETLib PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}") +target_link_libraries( + POETLib + PUBLIC RRuntime + PUBLIC IPhreeqcPOET + PUBLIC tug + PUBLIC MPI::MPI_C ) -target_compile_definitions(poetlib PUBLIC STRICT_R_HEADERS OMPI_SKIP_MPICXX) +# add_library(poetlib +# Base/Grid.cpp +# Base/SimParams.cpp +# Chemistry/ChemistryModule.cpp +# Chemistry/MasterFunctions.cpp +# Chemistry/WorkerFunctions.cpp +# Chemistry/SurrogateModels/DHT_Wrapper.cpp +# Chemistry/SurrogateModels/DHT.c +# Chemistry/SurrogateModels/HashFunctions.cpp +# Chemistry/SurrogateModels/InterpolationModule.cpp +# Chemistry/SurrogateModels/ProximityHashTable.cpp +# Transport/DiffusionModule.cpp +# ) -mark_as_advanced(PHREEQCRM_BUILD_MPI PHREEQCRM_DISABLE_OPENMP) -set(PHREEQCRM_DISABLE_OPENMP ON CACHE BOOL "" FORCE) +# target_link_libraries(poetlib PUBLIC +# DataStructures +# Init +# MPI::MPI_C +# ${MATH_LIBRARY} +# RRuntime +# PhreeqcRM +# tug +# ) -option(POET_DHT_DEBUG "Build with DHT debug info" OFF) +target_compile_definitions(POETLib PUBLIC STRICT_R_HEADERS OMPI_SKIP_MPICXX) -if(POET_DHT_DEBUG) - target_compile_definitions(poetlib PRIVATE DHT_STATISTICS) -endif() +# mark_as_advanced(PHREEQCRM_BUILD_MPI PHREEQCRM_DISABLE_OPENMP) +# set(PHREEQCRM_DISABLE_OPENMP ON CACHE BOOL "" FORCE) -option(POET_PHT_ADDITIONAL_INFO "Enables additional information in the PHT" OFF) +# option(POET_DHT_DEBUG "Build with DHT debug info" OFF) -if (POET_PHT_ADDITIONAL_INFO) - target_compile_definitions(poetlib PRIVATE POET_PHT_ADD) -endif() +# if(POET_DHT_DEBUG) +# target_compile_definitions(poetlib PRIVATE DHT_STATISTICS) +# endif() + +# option(POET_PHT_ADDITIONAL_INFO "Enables additional information in the PHT" OFF) + +# if (POET_PHT_ADDITIONAL_INFO) +# target_compile_definitions(poetlib PRIVATE POET_PHT_ADD) +# endif() file(READ "${PROJECT_SOURCE_DIR}/R_lib/kin_r_library.R" R_KIN_LIB ) +file(READ "${PROJECT_SOURCE_DIR}/R_lib/init_r_lib.R" R_INIT_LIB) configure_file(poet.hpp.in poet.hpp @ONLY) add_executable(poet poet.cpp) -target_link_libraries(poet PRIVATE poetlib MPI::MPI_C RRuntime) +target_link_libraries(poet PRIVATE POETLib MPI::MPI_C RRuntime) target_include_directories(poet PRIVATE "${CMAKE_CURRENT_BINARY_DIR}") -install(TARGETS poet DESTINATION bin) +add_executable(poet_init initializer.cpp) +target_link_libraries(poet_init PRIVATE POETLib RRuntime) +target_include_directories(poet_init PRIVATE "${CMAKE_CURRENT_BINARY_DIR}") + +install(TARGETS poet poet_init DESTINATION bin) + diff --git a/src/Chemistry/ChemistryDefs.hpp b/src/Chemistry/ChemistryDefs.hpp new file mode 100644 index 000000000..71c82edf6 --- /dev/null +++ b/src/Chemistry/ChemistryDefs.hpp @@ -0,0 +1,23 @@ +#pragma once + +#include +#include + +namespace poet { + +enum DHT_PROP_TYPES { DHT_TYPE_DEFAULT, DHT_TYPE_CHARGE, DHT_TYPE_TOTAL }; +enum CHEMISTRY_OUT_SOURCE { CHEM_PQC, CHEM_DHT, CHEM_INTERP }; + +struct WorkPackage { + std::size_t size; + std::vector> input; + std::vector> output; + std::vector mapping; + + WorkPackage(std::size_t _size) : size(_size) { + input.resize(size); + output.resize(size); + mapping.resize(size, CHEM_PQC); + } +}; +} // namespace poet \ No newline at end of file diff --git a/src/Chemistry/ChemistryModule.cpp b/src/Chemistry/ChemistryModule.cpp index 46531f2c6..155b90569 100644 --- a/src/Chemistry/ChemistryModule.cpp +++ b/src/Chemistry/ChemistryModule.cpp @@ -1,24 +1,17 @@ #include "ChemistryModule.hpp" +#include "PhreeqcEngine.hpp" #include "SurrogateModels/DHT_Wrapper.hpp" #include "SurrogateModels/Interpolation.hpp" -#include - #include #include -#include #include -#include #include #include -#include #include -#include #include -constexpr uint32_t MB_FACTOR = 1E6; - std::vector inverseDistanceWeighting(const std::vector &to_calc, const std::vector &from, @@ -26,20 +19,21 @@ inverseDistanceWeighting(const std::vector &to_calc, const std::vector> &output) { std::vector results = from; - const std::uint32_t buffer_size = input.size() + 1; - double buffer[buffer_size]; - double from_rescaled; + // const std::uint32_t buffer_size = input.size() + 1; + // double buffer[buffer_size]; + // double from_rescaled; const std::uint32_t data_set_n = input.size(); double rescaled[to_calc.size()][data_set_n + 1]; double weights[data_set_n]; // rescaling over all key elements - for (int key_comp_i = 0; key_comp_i < to_calc.size(); key_comp_i++) { + for (std::uint32_t key_comp_i = 0; key_comp_i < to_calc.size(); + key_comp_i++) { const auto output_comp_i = to_calc[key_comp_i]; // rescale input between 0 and 1 - for (int point_i = 0; point_i < data_set_n; point_i++) { + for (std::uint32_t point_i = 0; point_i < data_set_n; point_i++) { rescaled[key_comp_i][point_i] = input[point_i][key_comp_i]; } @@ -50,7 +44,7 @@ inverseDistanceWeighting(const std::vector &to_calc, const double max = *std::max_element(rescaled[key_comp_i], rescaled[key_comp_i] + data_set_n + 1); - for (int point_i = 0; point_i < data_set_n; point_i++) { + for (std::uint32_t point_i = 0; point_i < data_set_n; point_i++) { rescaled[key_comp_i][point_i] = ((max - min) != 0 ? (rescaled[key_comp_i][point_i] - min) / (max - min) @@ -62,9 +56,10 @@ inverseDistanceWeighting(const std::vector &to_calc, // calculate distances for each data set double inv_sum = 0; - for (int point_i = 0; point_i < data_set_n; point_i++) { + for (std::uint32_t point_i = 0; point_i < data_set_n; point_i++) { double distance = 0; - for (int key_comp_i = 0; key_comp_i < to_calc.size(); key_comp_i++) { + for (std::uint32_t key_comp_i = 0; key_comp_i < to_calc.size(); + key_comp_i++) { distance += std::pow( rescaled[key_comp_i][point_i] - rescaled[key_comp_i][data_set_n], 2); } @@ -79,7 +74,8 @@ inverseDistanceWeighting(const std::vector &to_calc, // bool has_h = false; // bool has_o = false; - for (int key_comp_i = 0; key_comp_i < to_calc.size(); key_comp_i++) { + for (std::uint32_t key_comp_i = 0; key_comp_i < to_calc.size(); + key_comp_i++) { const auto output_comp_i = to_calc[key_comp_i]; double key_delta = 0; @@ -91,7 +87,7 @@ inverseDistanceWeighting(const std::vector &to_calc, // has_o = true; // } - for (int j = 0; j < data_set_n; j++) { + for (std::uint32_t j = 0; j < data_set_n; j++) { key_delta += weights[j] * output[j][output_comp_i]; } @@ -158,25 +154,30 @@ inverseDistanceWeighting(const std::vector &to_calc, return results; } -poet::ChemistryModule::ChemistryModule(uint32_t nxyz, uint32_t wp_size, - std::uint32_t maxiter, - const ChemistryParams &chem_param, - MPI_Comm communicator) - : PhreeqcRM(nxyz, 1), group_comm(communicator), wp_size(wp_size), - params(chem_param) { +poet::ChemistryModule::ChemistryModule( + uint32_t wp_size_, const InitialList::ChemistryInit chem_params, + MPI_Comm communicator) + : group_comm(communicator), wp_size(wp_size_), params(chem_params) { + MPI_Comm_rank(communicator, &comm_rank); + MPI_Comm_size(communicator, &comm_size); - MPI_Comm_size(communicator, &this->comm_size); - MPI_Comm_rank(communicator, &this->comm_rank); + this->is_sequential = comm_size == 1; + this->is_master = comm_rank == 0; - this->is_sequential = (this->comm_size == 1); - this->is_master = (this->comm_rank == 0); + this->n_cells = chem_params.total_grid_cells; - if (!is_sequential && is_master) { - MPI_Bcast(&wp_size, 1, MPI_UINT32_T, 0, this->group_comm); - MPI_Bcast(&maxiter, 1, MPI_UINT32_T, 0, this->group_comm); + if (!is_master) { + for (const auto &[pqc_id, pqc_config] : chem_params.pqc_config) { + this->phreeqc_instances[pqc_id] = + std::make_unique(pqc_config); + } + // for (std::size_t i = 0; i < chem_params.pqc_ids.size(); i++) { + // this->phreeqc_instances[chem_params.pqc_ids[i]] = + // std::make_unique( + // chem_params.database, chem_params.pqc_scripts[i], + // chem_params.pqc_sol_order, chem_params.field_header, wp_size_); + // } } - - this->file_pad = std::ceil(std::log10(maxiter + 1)); } poet::ChemistryModule::~ChemistryModule() { @@ -185,168 +186,6 @@ poet::ChemistryModule::~ChemistryModule() { } } -poet::ChemistryModule -poet::ChemistryModule::createWorker(MPI_Comm communicator, - const ChemistryParams &chem_param) { - uint32_t wp_size; - MPI_Bcast(&wp_size, 1, MPI_UINT32_T, 0, communicator); - - std::uint32_t maxiter; - MPI_Bcast(&maxiter, 1, MPI_UINT32_T, 0, communicator); - - return ChemistryModule(wp_size, wp_size, maxiter, chem_param, communicator); -} - -void poet::ChemistryModule::RunInitFile(const std::string &input_script_path) { - if (this->is_master) { - int f_type = CHEM_INIT; - PropagateFunctionType(f_type); - - int count = input_script_path.size(); - ChemBCast(&count, 1, MPI_INT); - ChemBCast(const_cast(input_script_path.data()), count, MPI_CHAR); - } - - this->RunFile(true, true, false, input_script_path); - this->RunString(true, false, false, "DELETE; -all; PRINT; -warnings 0;"); - - this->FindComponents(); - - PhreeqcRM::initializePOET(this->speciesPerModule, this->prop_names); - this->prop_count = prop_names.size(); - - char exchange = (speciesPerModule[1] == 0 ? -1 : 1); - char kinetics = (speciesPerModule[2] == 0 ? -1 : 1); - char equilibrium = (speciesPerModule[3] == 0 ? -1 : 1); - char surface = (speciesPerModule[4] == 0 ? -1 : 1); - - std::vector ic1; - ic1.resize(this->nxyz * 7, -1); - // TODO: hardcoded reaction modules - for (int i = 0; i < nxyz; i++) { - ic1[i] = 1; // Solution 1 - ic1[nxyz + i] = equilibrium; // Equilibrium 1 - ic1[2 * nxyz + i] = exchange; // Exchange none - ic1[3 * nxyz + i] = surface; // Surface none - ic1[4 * nxyz + i] = -1; // Gas phase none - ic1[5 * nxyz + i] = -1; // Solid solutions none - ic1[6 * nxyz + i] = kinetics; // Kinetics 1 - } - - this->InitialPhreeqc2Module(ic1); -} - -void poet::ChemistryModule::initializeField(const Field &trans_field) { - - if (is_master) { - int f_type = CHEM_INIT_SPECIES; - PropagateFunctionType(f_type); - } - - std::vector essentials_backup{ - prop_names.begin() + speciesPerModule[0], prop_names.end()}; - - std::vector new_solution_names{ - this->prop_names.begin(), this->prop_names.begin() + speciesPerModule[0]}; - - if (is_master) { - for (auto &prop : trans_field.GetProps()) { - if (std::find(new_solution_names.begin(), new_solution_names.end(), - prop) == new_solution_names.end()) { - int size = prop.size(); - ChemBCast(&size, 1, MPI_INT); - ChemBCast(prop.data(), prop.size(), MPI_CHAR); - new_solution_names.push_back(prop); - } - } - int end = 0; - ChemBCast(&end, 1, MPI_INT); - } else { - constexpr int MAXSIZE = 128; - MPI_Status status; - int recv_size; - char recv_buffer[MAXSIZE]; - while (1) { - ChemBCast(&recv_size, 1, MPI_INT); - if (recv_size == 0) { - break; - } - ChemBCast(recv_buffer, recv_size, MPI_CHAR); - recv_buffer[recv_size] = '\0'; - new_solution_names.push_back(std::string(recv_buffer)); - } - } - - // now sort the new values - std::sort(new_solution_names.begin() + 3, new_solution_names.end()); - this->SetPOETSolutionNames(new_solution_names); - this->speciesPerModule[0] = new_solution_names.size(); - - // and append other processes than solutions - std::vector new_prop_names = new_solution_names; - new_prop_names.insert(new_prop_names.end(), essentials_backup.begin(), - essentials_backup.end()); - - std::vector old_prop_names{this->prop_names}; - this->prop_names = std::move(new_prop_names); - this->prop_count = prop_names.size(); - - if (is_master) { - this->n_cells = trans_field.GetRequestedVecSize(); - - std::vector> phreeqc_dump(this->nxyz); - this->getDumpedField(phreeqc_dump); - - if (is_sequential) { - std::vector init_vec; - for (std::size_t i = 0; i < n_cells; i++) { - init_vec.insert(init_vec.end(), phreeqc_dump[i].begin(), - phreeqc_dump[i].end()); - } - - const auto tmp_buffer{init_vec}; - this->unshuffleField(tmp_buffer, n_cells, prop_count, 1, init_vec); - this->chem_field = Field(n_cells, init_vec, prop_names); - return; - } - - std::vector &phreeqc_init = phreeqc_dump[0]; - std::vector> initial_values; - - for (const auto &vec : trans_field.As2DVector()) { - initial_values.push_back(vec); - } - - this->base_totals = {initial_values.at(0).at(0), - initial_values.at(1).at(0)}; - ChemBCast(base_totals.data(), 2, MPI_DOUBLE); - - for (int i = speciesPerModule[0]; i < phreeqc_init.size(); i++) { - std::vector init(n_cells, phreeqc_init[i]); - initial_values.push_back(std::move(init)); - } - - this->chem_field = Field(n_cells, initial_values, prop_names); - } else { - ChemBCast(base_totals.data(), 2, MPI_DOUBLE); - } - - if (this->params.use_dht || this->params.use_interp) { - initializeDHT(this->params.dht_size, this->params.dht_signifs); - setDHTSnapshots(this->params.dht_snaps, this->params.dht_outdir); - setDHTReadFile(this->params.dht_file); - - this->dht_enabled = this->params.use_dht; - - if (this->params.use_interp) { - initializeInterp(this->params.pht_max_entries, this->params.pht_size, - this->params.interp_min_entries, - this->params.pht_signifs); - this->interp_enabled = this->params.use_interp; - } - } -} - void poet::ChemistryModule::initializeDHT( uint32_t size_mb, const NamedVector &key_species) { constexpr uint32_t MB_FACTOR = 1E6; @@ -368,7 +207,7 @@ void poet::ChemistryModule::initializeDHT( if (key_species.empty()) { std::vector default_signif( - this->prop_names.size(), DHT_Wrapper::DHT_KEY_SIGNIF_DEFAULT); + this->prop_count, DHT_Wrapper::DHT_KEY_SIGNIF_DEFAULT); map_copy = NamedVector(this->prop_names, default_signif); } @@ -382,7 +221,7 @@ void poet::ChemistryModule::initializeDHT( this->dht = new DHT_Wrapper(dht_comm, dht_size, map_copy, key_indices, this->prop_names, params.hooks, - this->prop_count, params.use_interp); + this->prop_count, interp_enabled); this->dht->setBaseTotals(base_totals.at(0), base_totals.at(1)); } } @@ -393,6 +232,8 @@ inline std::vector poet::ChemistryModule::parseDHTSpeciesVec( std::vector species_indices; species_indices.reserve(key_species.size()); + const auto test = key_species.getNames(); + for (const auto &species : key_species.getNames()) { auto it = std::find(to_compare.begin(), to_compare.end(), species); if (it == to_compare.end()) { @@ -469,9 +310,9 @@ void poet::ChemistryModule::initializeInterp( if (key_species.empty()) { map_copy = this->dht->getKeySpecies(); - for (std::size_t i = 0; i < map_copy.size(); i++) { + for (auto i = 0; i < map_copy.size(); i++) { const std::uint32_t signif = - map_copy[i] - (map_copy[i] > InterpolationModule::COARSE_DIFF + static_cast(map_copy[i]) - (map_copy[i] > InterpolationModule::COARSE_DIFF ? InterpolationModule::COARSE_DIFF : 0); map_copy[i] = signif; diff --git a/src/Chemistry/ChemistryModule.hpp b/src/Chemistry/ChemistryModule.hpp index d343ffa31..7260422fb 100644 --- a/src/Chemistry/ChemistryModule.hpp +++ b/src/Chemistry/ChemistryModule.hpp @@ -2,16 +2,17 @@ #ifndef CHEMISTRYMODULE_H_ #define CHEMISTRYMODULE_H_ -#include "../Base/SimParams.hpp" -#include "../DataStructures/DataStructures.hpp" +#include "DataStructures/Field.hpp" +#include "DataStructures/NamedVector.hpp" + +#include "ChemistryDefs.hpp" + +#include "Init/InitialList.hpp" #include "SurrogateModels/DHT_Wrapper.hpp" #include "SurrogateModels/Interpolation.hpp" -#include -#include - +#include "PhreeqcEngine.hpp" #include -#include #include #include #include @@ -24,7 +25,7 @@ namespace poet { * \brief Wrapper around PhreeqcRM to provide POET specific parallelization with * easy access. */ -class ChemistryModule : public PhreeqcRM { +class ChemistryModule { public: /** * Creates a new instance of Chemistry module with given grid cell count, work @@ -41,61 +42,67 @@ public: * \param wp_size Count of grid cells to fill each work package at maximum. * \param communicator MPI communicator to distribute work in. */ - ChemistryModule(uint32_t nxyz, uint32_t wp_size, std::uint32_t maxiter, - const ChemistryParams &chem_param, MPI_Comm communicator); + ChemistryModule(uint32_t wp_size, + const InitialList::ChemistryInit chem_params, + MPI_Comm communicator); /** * Deconstructor, which frees DHT data structure if used. */ ~ChemistryModule(); - /** - * Parses the input script and extract information needed during runtime. - * - * **Only run by master**. - * - * Database must be loaded beforehand. - * - * \param input_script_path Path to input script to parse. - */ - void RunInitFile(const std::string &input_script_path); - + void masterSetField(Field field); /** * Run the chemical simulation with parameters set. */ - void RunCells(); + void simulate(double dt); - /** - * Returns the chemical field. - */ - auto GetField() const { return this->chem_field; } /** * Returns all known species names, including not only aqueous species, but * also equilibrium, exchange, surface and kinetic reactants. */ - auto GetPropNames() const { return this->prop_names; } + // auto GetPropNames() const { return this->prop_names; } /** * Return the accumulated runtime in seconds for chemical simulation. */ auto GetChemistryTime() const { return this->chem_t; } - /** - * Create a new worker instance inside given MPI communicator. - * - * Wraps communication needed before instanciation can take place. - * - * \param communicator MPI communicator to distribute work in. - * - * \returns A worker instance with fixed work package size. - */ - static ChemistryModule createWorker(MPI_Comm communicator, - const ChemistryParams &chem_param); + void setFilePadding(std::uint32_t maxiter) { + this->file_pad = + static_cast(std::ceil(std::log10(maxiter + 1))); + } - /** - * Default work package size. - */ - static constexpr uint32_t CHEM_DEFAULT_WP_SIZE = 5; + struct SurrogateSetup { + std::vector prop_names; + + bool dht_enabled; + std::uint32_t dht_size_mb; + + bool interp_enabled; + std::uint32_t interp_bucket_size; + std::uint32_t interp_size_mb; + std::uint32_t interp_min_entries; + }; + + void masterEnableSurrogates(const SurrogateSetup &setup) { + // FIXME: This is a hack to get the prop_names and prop_count from the setup + this->prop_names = setup.prop_names; + this->prop_count = setup.prop_names.size(); + + this->dht_enabled = setup.dht_enabled; + this->interp_enabled = setup.interp_enabled; + + if (this->dht_enabled || this->interp_enabled) { + this->initializeDHT(setup.dht_size_mb, this->params.dht_species); + } + + if (this->interp_enabled) { + this->initializeInterp(setup.interp_bucket_size, setup.interp_size_mb, + setup.interp_min_entries, + this->params.interp_species); + } + } /** * Intended to alias input parameters for grid initialization with a single @@ -118,17 +125,6 @@ public: DHT_SNAPS_ITEREND //!< output snapshots after each iteration }; - /** - * **This function has to be run!** - * - * Merge initial values from existing module with the chemistry module and set - * according internal variables. - * - * \param other Field to merge chemistry with. Most likely it is something - * like the diffusion field. - */ - void initializeField(const Field &other); - /** * **Only called by workers!** Start the worker listening loop. */ @@ -243,9 +239,7 @@ protected: const NamedVector &key_species); enum { - CHEM_INIT, - CHEM_WP_SIZE, - CHEM_INIT_SPECIES, + CHEM_FIELD_INIT, CHEM_DHT_ENABLE, CHEM_DHT_SIGNIF_VEC, CHEM_DHT_SNAPS, @@ -294,7 +288,7 @@ protected: using worker_list_t = std::vector; using workpointer_t = std::vector::iterator; - void MasterRunParallel(); + void MasterRunParallel(double dt); void MasterRunSequential(); void MasterSendPkgs(worker_list_t &w_list, workpointer_t &work_pointer, @@ -320,8 +314,8 @@ protected: void WorkerPerfToMaster(int type, const struct worker_s &timings); void WorkerMetricsToMaster(int type); - IRM_RESULT WorkerRunWorkPackage(WorkPackage &work_package, double dSimTime, - double dTimestep); + void WorkerRunWorkPackage(WorkPackage &work_package, double dSimTime, + double dTimestep); std::vector CalculateWPSizesVector(uint32_t n_cells, uint32_t wp_size) const; @@ -344,7 +338,7 @@ protected: bool is_sequential; bool is_master; - uint32_t wp_size{CHEM_DEFAULT_WP_SIZE}; + uint32_t wp_size; bool dht_enabled{false}; int dht_snaps_type{DHT_SNAPS_DISABLED}; std::string dht_file_out_dir; @@ -372,7 +366,7 @@ protected: bool print_progessbar{false}; - std::uint32_t file_pad; + std::uint8_t file_pad{1}; double chem_t{0.}; @@ -382,11 +376,9 @@ protected: Field chem_field; - static constexpr int MODULE_COUNT = 5; + const InitialList::ChemistryInit params; - const ChemistryParams ¶ms; - - std::array speciesPerModule{}; + std::map> phreeqc_instances; }; } // namespace poet diff --git a/src/Chemistry/MasterFunctions.cpp b/src/Chemistry/MasterFunctions.cpp index 543cd21a8..9ee52be82 100644 --- a/src/Chemistry/MasterFunctions.cpp +++ b/src/Chemistry/MasterFunctions.cpp @@ -1,12 +1,9 @@ #include "ChemistryModule.hpp" -#include - #include #include #include #include -#include #include std::vector @@ -308,45 +305,46 @@ inline void poet::ChemistryModule::MasterRecvPkgs(worker_list_t &w_list, } } -void poet::ChemistryModule::RunCells() { +void poet::ChemistryModule::simulate(double dt) { double start_t{MPI_Wtime()}; if (this->is_sequential) { MasterRunSequential(); return; } - MasterRunParallel(); + MasterRunParallel(dt); double end_t{MPI_Wtime()}; this->chem_t += end_t - start_t; } void poet::ChemistryModule::MasterRunSequential() { - std::vector shuffled_field = - shuffleField(chem_field.AsVector(), n_cells, prop_count, 1); + // std::vector shuffled_field = + // shuffleField(chem_field.AsVector(), n_cells, prop_count, 1); - std::vector> input; - for (std::size_t i = 0; i < n_cells; i++) { - input.push_back( - std::vector(shuffled_field.begin() + (i * prop_count), - shuffled_field.begin() + ((i + 1) * prop_count))); - } + // std::vector> input; + // for (std::size_t i = 0; i < n_cells; i++) { + // input.push_back( + // std::vector(shuffled_field.begin() + (i * prop_count), + // shuffled_field.begin() + ((i + 1) * + // prop_count))); + // } - this->setDumpedField(input); - PhreeqcRM::RunCells(); - this->getDumpedField(input); + // this->setDumpedField(input); + // PhreeqcRM::RunCells(); + // this->getDumpedField(input); - shuffled_field.clear(); - for (std::size_t i = 0; i < n_cells; i++) { - shuffled_field.insert(shuffled_field.end(), input[i].begin(), - input[i].end()); - } + // shuffled_field.clear(); + // for (std::size_t i = 0; i < n_cells; i++) { + // shuffled_field.insert(shuffled_field.end(), input[i].begin(), + // input[i].end()); + // } - std::vector out_vec{shuffled_field}; - unshuffleField(shuffled_field, n_cells, prop_count, 1, out_vec); - chem_field = out_vec; + // std::vector out_vec{shuffled_field}; + // unshuffleField(shuffled_field, n_cells, prop_count, 1, out_vec); + // chem_field = out_vec; } -void poet::ChemistryModule::MasterRunParallel() { +void poet::ChemistryModule::MasterRunParallel(double dt) { /* declare most of the needed variables here */ double seq_a, seq_b, seq_c, seq_d; double worker_chemistry_a, worker_chemistry_b; @@ -360,7 +358,6 @@ void poet::ChemistryModule::MasterRunParallel() { MPI_Barrier(this->group_comm); - double dt = this->PhreeqcRM::GetTimeStep(); static uint32_t iteration = 0; /* start time measurement of sequential part */ @@ -466,3 +463,13 @@ poet::ChemistryModule::CalculateWPSizesVector(uint32_t n_cells, return wp_sizes_vector; } + +void poet::ChemistryModule::masterSetField(Field field) { + this->chem_field = field; + this->prop_count = field.GetProps().size(); + + int ftype = CHEM_FIELD_INIT; + PropagateFunctionType(ftype); + + ChemBCast(&this->prop_count, 1, MPI_UINT32_T); +} \ No newline at end of file diff --git a/src/Chemistry/SurrogateModels/DHT_Wrapper.cpp b/src/Chemistry/SurrogateModels/DHT_Wrapper.cpp index 69491b2c2..eac73b14a 100644 --- a/src/Chemistry/SurrogateModels/DHT_Wrapper.cpp +++ b/src/Chemistry/SurrogateModels/DHT_Wrapper.cpp @@ -22,13 +22,15 @@ #include "DHT_Wrapper.hpp" +#include "Init/InitialList.hpp" +#include "Rounding.hpp" + #include #include #include #include #include #include -#include #include #include @@ -40,7 +42,7 @@ DHT_Wrapper::DHT_Wrapper(MPI_Comm dht_comm, std::uint64_t dht_size, const NamedVector &key_species, const std::vector &key_indices, const std::vector &_output_names, - const ChemistryParams::Chem_Hook_Functions &_hooks, + const InitialList::ChemistryHookFunctions &_hooks, uint32_t data_count, bool _with_interp) : key_count(key_indices.size()), data_count(data_count), input_key_elements(key_indices), communicator(dht_comm), diff --git a/src/Chemistry/SurrogateModels/DHT_Wrapper.hpp b/src/Chemistry/SurrogateModels/DHT_Wrapper.hpp index 791be2d1b..fed271923 100644 --- a/src/Chemistry/SurrogateModels/DHT_Wrapper.hpp +++ b/src/Chemistry/SurrogateModels/DHT_Wrapper.hpp @@ -23,19 +23,18 @@ #ifndef DHT_WRAPPER_H #define DHT_WRAPPER_H -#include "../../Base/RInsidePOET.hpp" -#include "../../Base/SimParams.hpp" -#include "../../DataStructures/DataStructures.hpp" -#include "../enums.hpp" -#include "HashFunctions.hpp" +#include "Base/RInsidePOET.hpp" +#include "DataStructures/NamedVector.hpp" + +#include "Chemistry/ChemistryDefs.hpp" + +#include "Init/InitialList.hpp" #include "LookupKey.hpp" -#include "Rounding.hpp" #include #include #include #include -#include #include #include @@ -88,7 +87,7 @@ public: const NamedVector &key_species, const std::vector &key_indices, const std::vector &output_names, - const ChemistryParams::Chem_Hook_Functions &hooks, + const InitialList::ChemistryHookFunctions &hooks, uint32_t data_count, bool with_interp); /** * @brief Destroy the dht wrapper object @@ -260,7 +259,7 @@ private: const std::vector &output_names; - const ChemistryParams::Chem_Hook_Functions &hooks; + const InitialList::ChemistryHookFunctions &hooks; const bool with_interp; DHT_ResultObject dht_results; diff --git a/src/Chemistry/SurrogateModels/Interpolation.hpp b/src/Chemistry/SurrogateModels/Interpolation.hpp index 15361f4bf..a1886965d 100644 --- a/src/Chemistry/SurrogateModels/Interpolation.hpp +++ b/src/Chemistry/SurrogateModels/Interpolation.hpp @@ -3,14 +3,13 @@ #ifndef INTERPOLATION_H_ #define INTERPOLATION_H_ -#include "../../Base/SimParams.hpp" -#include "../../DataStructures/DataStructures.hpp" +#include "DataStructures/NamedVector.hpp" + #include "DHT_Wrapper.hpp" +#include "Init/InitialList.hpp" #include "LookupKey.hpp" #include "Rounding.hpp" -#include -#include #include #include #include @@ -22,7 +21,6 @@ extern "C" { } #include -#include #include #include @@ -165,7 +163,7 @@ public: const NamedVector &interp_key_signifs, const std::vector &dht_key_indices, const std::vector &out_names, - const ChemistryParams::Chem_Hook_Functions &hooks); + const InitialList::ChemistryHookFunctions &hooks); enum result_status { RES_OK, INSUFFICIENT_DATA, NOT_NEEDED }; @@ -261,7 +259,7 @@ private: return out_key; } - const ChemistryParams::Chem_Hook_Functions &hooks; + const InitialList::ChemistryHookFunctions &hooks; const std::vector &out_names; const std::vector dht_names; }; diff --git a/src/Chemistry/SurrogateModels/InterpolationModule.cpp b/src/Chemistry/SurrogateModels/InterpolationModule.cpp index 26bfd4ab0..455c96729 100644 --- a/src/Chemistry/SurrogateModels/InterpolationModule.cpp +++ b/src/Chemistry/SurrogateModels/InterpolationModule.cpp @@ -1,8 +1,9 @@ // Time-stamp: "Last modified 2023-08-16 17:02:31 mluebke" +#include "Init/InitialList.hpp" #include "Interpolation.hpp" -#include "../../DataStructures/DataStructures.hpp" #include "DHT_Wrapper.hpp" +#include "DataStructures/NamedVector.hpp" #include "HashFunctions.hpp" #include "LookupKey.hpp" #include "Rounding.hpp" @@ -35,7 +36,7 @@ InterpolationModule::InterpolationModule( const NamedVector &interp_key_signifs, const std::vector &dht_key_indices, const std::vector &_out_names, - const ChemistryParams::Chem_Hook_Functions &_hooks) + const InitialList::ChemistryHookFunctions &_hooks) : dht_instance(dht), key_signifs(interp_key_signifs), key_indices(dht_key_indices), min_entries_needed(min_entries_needed), dht_names(dht.getKeySpecies().getNames()), out_names(_out_names), diff --git a/src/Chemistry/WorkerFunctions.cpp b/src/Chemistry/WorkerFunctions.cpp index 4e67e33d5..96a67898d 100644 --- a/src/Chemistry/WorkerFunctions.cpp +++ b/src/Chemistry/WorkerFunctions.cpp @@ -2,10 +2,9 @@ #include "SurrogateModels/DHT_Wrapper.hpp" #include "SurrogateModels/Interpolation.hpp" -#include +#include "Chemistry/ChemistryDefs.hpp" + #include -#include -#include #include #include #include @@ -17,6 +16,7 @@ #include namespace poet { + inline std::string get_string(int root, MPI_Comm communicator) { int count; MPI_Bcast(&count, 1, MPI_INT, root, communicator); @@ -45,13 +45,8 @@ void poet::ChemistryModule::WorkerLoop() { PropagateFunctionType(func_type); switch (func_type) { - case CHEM_INIT: { - RunInitFile(get_string(0, this->group_comm)); - break; - } - case CHEM_INIT_SPECIES: { - Field dummy; - initializeField(dummy); + case CHEM_FIELD_INIT: { + ChemBCast(&this->prop_count, 1, MPI_UINT32_T); break; } case CHEM_WORK_LOOP: { @@ -177,9 +172,7 @@ void poet::ChemistryModule::WorkerDoWork(MPI_Status &probe_status, phreeqc_time_start = MPI_Wtime(); - if (WorkerRunWorkPackage(s_curr_wp, current_sim_time, dt) != IRM_OK) { - std::cerr << "Phreeqc error" << std::endl; - }; + WorkerRunWorkPackage(s_curr_wp, current_sim_time, dt); phreeqc_time_end = MPI_Wtime(); @@ -285,42 +278,34 @@ void poet::ChemistryModule::WorkerReadDHTDump( } } -IRM_RESULT -poet::ChemistryModule::WorkerRunWorkPackage(WorkPackage &work_package, - double dSimTime, double dTimestep) { - // check if we actually need to start phreeqc - std::vector pqc_mapping; +void poet::ChemistryModule::WorkerRunWorkPackage(WorkPackage &work_package, + double dSimTime, + double dTimestep) { - for (std::size_t i = 0; i < work_package.size; i++) { - if (work_package.mapping[i] == CHEM_PQC) { - pqc_mapping.push_back(i); + for (std::size_t wp_id = 0; wp_id < work_package.size; wp_id++) { + if (work_package.mapping[wp_id] != CHEM_PQC) { + continue; + } + + auto curr_input = work_package.input[wp_id]; + const auto pqc_id = static_cast(curr_input[0]); + + auto &phreeqc_instance = this->phreeqc_instances[pqc_id]; + work_package.output[wp_id] = work_package.input[wp_id]; + + curr_input.erase(std::remove_if(curr_input.begin(), curr_input.end(), + [](double d) { return std::isnan(d); }), + curr_input.end()); + + phreeqc_instance->runCell(curr_input, dTimestep); + + std::size_t output_index = 0; + for (std::size_t i = 0; i < work_package.output[wp_id].size(); i++) { + if (std::isnan(work_package.output[wp_id][i])) { + work_package.output[wp_id][i] = curr_input[output_index++]; + } } } - - if (pqc_mapping.empty()) { - return IRM_OK; - } - - IRM_RESULT result; - this->PhreeqcRM::setPOETMapping(pqc_mapping); - this->setDumpedField(work_package.input); - - this->PhreeqcRM::SetTime(dSimTime); - this->PhreeqcRM::SetTimeStep(dTimestep); - - result = this->PhreeqcRM::RunCells(); - - std::vector> output_tmp(work_package.size); - - this->getDumpedField(output_tmp); - - for (std::size_t i = 0; i < work_package.size; i++) { - if (work_package.mapping[i] == CHEM_PQC) { - work_package.output[i] = output_tmp[i]; - } - } - - return result; } void poet::ChemistryModule::WorkerPerfToMaster(int type, diff --git a/src/Chemistry/enums.hpp b/src/Chemistry/enums.hpp deleted file mode 100644 index 11c3c1281..000000000 --- a/src/Chemistry/enums.hpp +++ /dev/null @@ -1,10 +0,0 @@ -#ifndef ENUMS_H_ -#define ENUMS_H_ - -namespace poet { -enum DHT_PROP_TYPES { DHT_TYPE_DEFAULT, DHT_TYPE_CHARGE, DHT_TYPE_TOTAL }; - -enum CHEMISTRY_OUT_SOURCE { CHEM_PQC, CHEM_DHT, CHEM_INTERP }; -} // namespace poet - -#endif // ENUMS_H_ diff --git a/src/DataStructures/CMakeLists.txt b/src/DataStructures/CMakeLists.txt new file mode 100644 index 000000000..bbc0263a8 --- /dev/null +++ b/src/DataStructures/CMakeLists.txt @@ -0,0 +1,3 @@ +target_sources(POETLib PRIVATE + Field.cpp +) \ No newline at end of file diff --git a/src/DataStructures/Field.cpp b/src/DataStructures/Field.cpp index 71d51dea4..6eb6cab1f 100644 --- a/src/DataStructures/Field.cpp +++ b/src/DataStructures/Field.cpp @@ -1,6 +1,9 @@ -#include "DataStructures.hpp" +#include "Field.hpp" #include +#include +#include +#include #include #include #include @@ -117,40 +120,14 @@ poet::FieldColumn &poet::Field::operator[](const std::string &key) { } SEXP poet::Field::asSEXP() const { - const std::size_t cols = this->props.size(); + Rcpp::List output; - SEXP s_names = PROTECT(Rf_allocVector(STRSXP, cols)); - SEXP s_output = PROTECT(Rf_allocVector(VECSXP, cols)); - - for (std::size_t prop_i = 0; prop_i < this->props.size(); prop_i++) { - const auto &name = this->props[prop_i]; - SEXP s_values = PROTECT(Rf_allocVector(REALSXP, this->req_vec_size)); - const auto values = this->find(name)->second; - - SET_STRING_ELT(s_names, prop_i, Rf_mkChar(name.c_str())); - - for (std::size_t i = 0; i < this->req_vec_size; i++) { - REAL(s_values)[i] = values[i]; - } - - SET_VECTOR_ELT(s_output, prop_i, s_values); - - UNPROTECT(1); + for (const auto &elem : this->props) { + const auto map_it = this->find(elem); + output[elem] = Rcpp::wrap(map_it->second); } - SEXP s_rownames = PROTECT(Rf_allocVector(INTSXP, this->req_vec_size)); - for (std::size_t i = 0; i < this->req_vec_size; i++) { - INTEGER(s_rownames)[i] = static_cast(i + 1); - } - - Rf_setAttrib(s_output, R_ClassSymbol, - Rf_ScalarString(Rf_mkChar("data.frame"))); - Rf_setAttrib(s_output, R_NamesSymbol, s_names); - Rf_setAttrib(s_output, R_RowNamesSymbol, s_rownames); - - UNPROTECT(3); - - return s_output; + return output; } poet::Field &poet::Field::operator=(const FieldColumn &cont_field) { @@ -199,34 +176,26 @@ poet::Field::operator=(const std::vector &cont_field) { void poet::Field::fromSEXP(const SEXP &s_rhs) { this->clear(); - SEXP s_vec = PROTECT(Rf_coerceVector(s_rhs, VECSXP)); - SEXP s_names = PROTECT(Rf_getAttrib(s_vec, R_NamesSymbol)); + Rcpp::List in_list; - std::size_t cols = static_cast(Rf_length(s_vec)); - - this->props.clear(); - this->props.reserve(cols); - - for (std::size_t i = 0; i < cols; i++) { - const std::string prop_name(CHAR(STRING_ELT(s_names, i))); - this->props.push_back(prop_name); - - SEXP s_values = PROTECT(VECTOR_ELT(s_vec, i)); - - if (i == 0) { - this->req_vec_size = static_cast(Rf_length(s_values)); - } - - FieldColumn input(this->req_vec_size); - - for (std::size_t j = 0; j < this->req_vec_size; j++) { - input[j] = static_cast(REAL(s_values)[j]); - } - - UNPROTECT(1); - - this->insert({prop_name, input}); + try { + in_list = Rcpp::List(s_rhs); + } catch (const Rcpp::exception &e) { + throw std::runtime_error("Input cannot be casted as list."); } - UNPROTECT(2); + if (in_list.size() == 0) { + return; + } + + this->props = Rcpp::as>(in_list.names()); + + const Rcpp::NumericVector &in_vec = in_list[0]; + + this->req_vec_size = static_cast(in_vec.size()); + + for (const auto &elem : this->props) { + const auto values = Rcpp::as>(in_list[elem]); + this->insert({elem, values}); + } } diff --git a/src/DataStructures/DataStructures.hpp b/src/DataStructures/Field.hpp similarity index 84% rename from src/DataStructures/DataStructures.hpp rename to src/DataStructures/Field.hpp index 30f69e5dd..90119f2a6 100644 --- a/src/DataStructures/DataStructures.hpp +++ b/src/DataStructures/Field.hpp @@ -1,54 +1,15 @@ #ifndef DATASTRUCTURES_H_ #define DATASTRUCTURES_H_ -#include "../Chemistry/enums.hpp" - #include -#include -#include #include -#include -#include +#include #include -#include #include namespace poet { -struct WorkPackage { - std::size_t size; - std::vector> input; - std::vector> output; - std::vector mapping; - - WorkPackage(size_t _size) : size(_size) { - input.resize(size); - output.resize(size); - mapping.resize(size, CHEM_PQC); - } -}; - -template class NamedVector : public Rcpp::NumericVector { -public: - NamedVector() : Rcpp::NumericVector(){}; - - NamedVector(const std::vector &in_names, - const std::vector &in_values) - : Rcpp::NumericVector(Rcpp::wrap(in_values)) { - this->names() = Rcpp::CharacterVector(Rcpp::wrap(in_names)); - } - - NamedVector(const SEXP &s) : Rcpp::NumericVector(s){}; - - bool empty() const { return (this->size() == 0); } - - std::vector getNames() const { - return Rcpp::as>(this->names()); - } - std::vector getValues() const { return Rcpp::as>(*this); } -}; - using FieldColumn = std::vector; /** diff --git a/src/DataStructures/NamedVector.hpp b/src/DataStructures/NamedVector.hpp new file mode 100644 index 000000000..f120b0e48 --- /dev/null +++ b/src/DataStructures/NamedVector.hpp @@ -0,0 +1,25 @@ +#pragma once + +#include + +namespace poet { +template class NamedVector : public Rcpp::NumericVector { +public: + NamedVector() : Rcpp::NumericVector(){}; + + NamedVector(const std::vector &in_names, + const std::vector &in_values) + : Rcpp::NumericVector(Rcpp::wrap(in_values)) { + this->names() = Rcpp::CharacterVector(Rcpp::wrap(in_names)); + } + + NamedVector(const SEXP &s) : Rcpp::NumericVector(s){}; + + bool empty() const { return (this->size() == 0); } + + std::vector getNames() const { + return Rcpp::as>(this->names()); + } + std::vector getValues() const { return Rcpp::as>(*this); } +}; +} // namespace poet \ No newline at end of file diff --git a/src/Init/ChemistryInit.cpp b/src/Init/ChemistryInit.cpp new file mode 100644 index 000000000..5eb0d01f4 --- /dev/null +++ b/src/Init/ChemistryInit.cpp @@ -0,0 +1,86 @@ +#include "InitialList.hpp" + +#include +#include +#include + +namespace poet { + +void InitialList::initChemistry(const Rcpp::List &chem) { + this->pqc_solutions = std::vector( + this->transport_names.begin() + 3, this->transport_names.end()); + + this->pqc_solution_primaries = this->phreeqc->getSolutionPrimaries(); + + if (chem.containsElementNamed("dht_species")) { + this->dht_species = Rcpp::as>(chem["dht_species"]); + } + + if (chem.containsElementNamed("pht_species")) { + this->interp_species = Rcpp::as>(chem["pht_species"]); + } + + if (chem.containsElementNamed("hooks")) { + this->chem_hooks = Rcpp::List(chem["hooks"]); + + std::vector hook_names = this->chem_hooks.names(); + + for (const auto &name : hook_names) { + if (this->hook_name_list.find(name) == this->hook_name_list.end()) { + Rcpp::Rcerr << "Unknown chemistry hook: " << name << std::endl; + } + } + } + + this->field_header = + Rcpp::as>(this->initial_grid.names()); + this->field_header.erase(this->field_header.begin()); +} + +InitialList::ChemistryInit InitialList::getChemistryInit() const { + ChemistryInit chem_init; + + chem_init.total_grid_cells = this->n_cols * this->n_rows; + + // chem_init.field_header = this->field_header; + + chem_init.database = database; + // chem_init.pqc_scripts = pqc_scripts; + // chem_init.pqc_ids = pqc_ids; + + for (std::size_t i = 0; i < pqc_scripts.size(); i++) { + POETInitCell cell = { + pqc_solutions, + pqc_solution_primaries, + Rcpp::as>(pqc_exchanger[i]), + Rcpp::as>(pqc_kinetics[i]), + Rcpp::as>(pqc_equilibrium[i]), + Rcpp::as>(pqc_surface_comps[i]), + Rcpp::as>(pqc_surface_charges[i])}; + + chem_init.pqc_config[pqc_ids[i]] = {database, pqc_scripts[i], cell}; + } + + // chem_init.pqc_sol_order = pqc_solutions; + + chem_init.dht_species = dht_species; + chem_init.interp_species = interp_species; + + if (this->chem_hooks.size() > 0) { + if (this->chem_hooks.containsElementNamed("dht_fill")) { + chem_init.hooks.dht_fill = this->chem_hooks["dht_fill"]; + } + if (this->chem_hooks.containsElementNamed("dht_fuzz")) { + chem_init.hooks.dht_fuzz = this->chem_hooks["dht_fuzz"]; + } + if (this->chem_hooks.containsElementNamed("interp_pre")) { + chem_init.hooks.interp_pre = this->chem_hooks["interp_pre"]; + } + if (this->chem_hooks.containsElementNamed("interp_post")) { + chem_init.hooks.interp_post = this->chem_hooks["interp_post"]; + } + } + + return chem_init; +} +} // namespace poet \ No newline at end of file diff --git a/src/Init/DiffusionInit.cpp b/src/Init/DiffusionInit.cpp new file mode 100644 index 000000000..3eb2fe2fd --- /dev/null +++ b/src/Init/DiffusionInit.cpp @@ -0,0 +1,406 @@ +#include +#include +// leave above Rcpp includes, as eigen seem to have problems with a preceding +// Rcpp include + +#include +#include + +#include "DataStructures/Field.hpp" +#include "InitialList.hpp" +#include "PhreeqcInit.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace poet { + +enum SEXP_TYPE { SEXP_IS_LIST, SEXP_IS_VEC }; + +const std::map tug_side_mapping = { + {tug::BC_SIDE_RIGHT, "E"}, + {tug::BC_SIDE_LEFT, "W"}, + {tug::BC_SIDE_TOP, "N"}, + {tug::BC_SIDE_BOTTOM, "S"}}; + +static std::vector colMajToRowMaj(const Rcpp::NumericVector &vec, + std::uint32_t n_cols, + std::uint32_t n_rows) { + if (vec.size() == 1) { + return std::vector(n_cols * n_rows, vec[0]); + } else { + if (vec.size() != n_cols * n_rows) { + throw std::runtime_error("Alpha matrix does not match grid dimensions"); + } + + std::vector alpha(n_cols * n_rows); + + for (std::uint32_t i = 0; i < n_cols; i++) { + for (std::uint32_t j = 0; j < n_rows; j++) { + alpha[i * n_rows + j] = vec[j * n_cols + i]; + } + } + return alpha; + } +} + +static std::vector +extend_transport_names(std::unique_ptr &phreeqc, + const Rcpp::List &boundaries_list, + const Rcpp::List &inner_boundaries, + const std::vector &old_trans_names) { + + std::vector transport_names = old_trans_names; + std::set constant_pqc_ids; + + for (const auto &side : tug_side_mapping) { + if (!boundaries_list.containsElementNamed(side.second.c_str())) { + continue; + } + + Rcpp::List mapping = boundaries_list[side.second]; + + const Rcpp::NumericVector cells = mapping["cell"]; + const Rcpp::NumericVector values = mapping["sol_id"]; + const Rcpp::CharacterVector type_str = mapping["type"]; + + if (cells.size() != values.size()) { + throw std::runtime_error("Boundary [" + side.second + + "] cells and values are not the same " + "length"); + } + + for (auto i = 0; i < cells.size(); i++) { + if (type_str[i] == "constant") { + constant_pqc_ids.insert(static_cast(values[i])); + } + } + } + + if (inner_boundaries.size() > 0) { + const Rcpp::NumericVector values = inner_boundaries["sol_id"]; + for (auto i = 0; i < values.size(); i++) { + constant_pqc_ids.insert(static_cast(values[i])); + } + } + + if (!constant_pqc_ids.empty()) { + constexpr std::size_t keep_h_o_charge = 3; + + for (const auto &pqc_id : constant_pqc_ids) { + const auto solution_names = phreeqc->getSolutionNames(pqc_id); + + // add those strings which are not already in the transport_names + for (const auto &name : solution_names) { + if (std::find(transport_names.begin(), transport_names.end(), name) == + transport_names.end()) { + auto position = + std::lower_bound(transport_names.begin() + keep_h_o_charge, + transport_names.end(), name); + + transport_names.insert(position, name); + } + } + } + } + + return transport_names; +} + +static Rcpp::List +extend_initial_grid(const Rcpp::List &initial_grid, + const std::vector &transport_names) { + // std::vector names_to_add(transport_names.begin() + old_size, + // transport_names.end()); + + Rcpp::Function extend_grid_R("add_missing_transport_species"); + + return extend_grid_R(initial_grid, Rcpp::wrap(transport_names)); +} + +std::pair +InitialList::resolveBoundaries(const Rcpp::List &boundaries_list, + const Rcpp::List &inner_boundaries) { + Rcpp::List bound_list; + Rcpp::List inner_bound; + Rcpp::Function resolve_R("resolve_pqc_bound"); + + const std::size_t old_transport_size = this->transport_names.size(); + + this->transport_names = extend_transport_names( + this->phreeqc, boundaries_list, inner_boundaries, this->transport_names); + + const std::size_t new_transport_size = this->transport_names.size(); + + if (old_transport_size != new_transport_size) { + this->initial_grid = + extend_initial_grid(this->initial_grid, this->transport_names); + } + + for (const auto &species : this->transport_names) { + Rcpp::List spec_list; + + for (const auto &side : tug_side_mapping) { + Rcpp::NumericVector c_type = + (side.first == tug::BC_SIDE_RIGHT || side.first == tug::BC_SIDE_LEFT) + ? Rcpp::NumericVector(n_rows, tug::BC_TYPE_CLOSED) + : Rcpp::NumericVector(n_cols, tug::BC_TYPE_CLOSED); + + Rcpp::NumericVector c_value = + (side.first == tug::BC_SIDE_RIGHT || side.first == tug::BC_SIDE_LEFT) + ? Rcpp::NumericVector(n_rows, 0) + : Rcpp::NumericVector(n_cols, 0); + + if (boundaries_list.containsElementNamed(side.second.c_str())) { + const Rcpp::List mapping = boundaries_list[side.second]; + + const Rcpp::NumericVector cells = mapping["cell"]; + const Rcpp::NumericVector values = mapping["sol_id"]; + const Rcpp::CharacterVector type_str = mapping["type"]; + + if (cells.size() != values.size()) { + throw std::runtime_error( + "Boundary [" + side.second + + "] cells and values are not the same length"); + } + + for (auto i = 0; i < cells.size(); i++) { + const auto c_id = cells[i] - 1; + if (type_str[i] == "closed") { + c_type[c_id] = tug::BC_TYPE_CLOSED; + } else if (type_str[i] == "constant") { + c_type[c_id] = tug::BC_TYPE_CONSTANT; + c_value[c_id] = Rcpp::as( + resolve_R(this->phreeqc_mat, Rcpp::wrap(species), values[i])); + } else { + throw std::runtime_error("Unknown boundary type"); + } + } + } + + spec_list[side.second] = Rcpp::List::create( + Rcpp::Named("type") = c_type, Rcpp::Named("value") = c_value); + } + + bound_list[species] = spec_list; + + if (inner_boundaries.size() > 0) { + const Rcpp::NumericVector &inner_row_vec = inner_boundaries["row"]; + const Rcpp::NumericVector &inner_col_vec = inner_boundaries["col"]; + const Rcpp::NumericVector &inner_pqc_id_vec = inner_boundaries["sol_id"]; + + std::vector rows; + std::vector cols; + std::vector c_value; + + if (inner_row_vec.size() != inner_col_vec.size() || + inner_row_vec.size() != inner_pqc_id_vec.size()) { + throw std::runtime_error( + "Inner boundary vectors are not the same length"); + } + + for (std::size_t i = 0; i < inner_row_vec.size(); i++) { + rows.push_back(inner_row_vec[i] - 1); + cols.push_back(inner_col_vec[i] - 1); + c_value.push_back(Rcpp::as(resolve_R( + this->phreeqc_mat, Rcpp::wrap(species), inner_pqc_id_vec[i]))); + } + + inner_bound[species] = + Rcpp::List::create(Rcpp::Named("row") = Rcpp::wrap(rows), + Rcpp::Named("col") = Rcpp::wrap(cols), + Rcpp::Named("value") = Rcpp::wrap(c_value)); + } + } + + return std::make_pair(bound_list, inner_bound); +} + +static inline SEXP_TYPE get_datatype(const SEXP &input) { + Rcpp::Function check_list("is.list"); + + if (Rcpp::as(check_list(input))) { + return SEXP_IS_LIST; + } else { + return SEXP_IS_VEC; + } +} + +static Rcpp::List parseAlphas(const SEXP &input, + const std::vector &transport_names, + std::uint32_t n_cols, std::uint32_t n_rows) { + Rcpp::List out_list; + + SEXP_TYPE input_type = get_datatype(input); + + switch (input_type) { + case SEXP_IS_LIST: { + const Rcpp::List input_list(input); + + if (input_list.size() != transport_names.size()) { + throw std::runtime_error( + "Length of alphas list does not match number of transport names"); + } + + for (const auto &name : transport_names) { + if (!input_list.containsElementNamed(name.c_str())) { + throw std::runtime_error("Alphas list does not contain transport name"); + } + + const Rcpp::NumericVector &alpha_col_order_vec = input_list[name]; + + out_list[name] = + Rcpp::wrap(colMajToRowMaj(alpha_col_order_vec, n_cols, n_rows)); + } + break; + } + case SEXP_IS_VEC: { + const Rcpp::NumericVector alpha(input); + for (const auto &name : transport_names) { + out_list[name] = Rcpp::wrap(colMajToRowMaj(alpha, n_cols, n_rows)); + } + break; + } + default: + throw std::runtime_error("Unknown data type"); + } + + return out_list; +} +void InitialList::initDiffusion(const Rcpp::List &diffusion_input) { + Rcpp::List boundaries; + Rcpp::List inner_boundaries; + + if (diffusion_input.containsElementNamed( + DIFFU_MEMBER_STR(DiffusionMembers::BOUNDARIES))) { + boundaries = + diffusion_input[DIFFU_MEMBER_STR(DiffusionMembers::BOUNDARIES)]; + } + + if (diffusion_input.containsElementNamed( + DIFFU_MEMBER_STR(DiffusionMembers::INNER_BOUNDARIES))) { + inner_boundaries = + diffusion_input[DIFFU_MEMBER_STR(DiffusionMembers::INNER_BOUNDARIES)]; + } + + const SEXP &alpha_x = + diffusion_input[DIFFU_MEMBER_STR(DiffusionMembers::ALPHA_X)]; + const SEXP &alpha_y = + diffusion_input[DIFFU_MEMBER_STR(DiffusionMembers::ALPHA_Y)]; + + const auto resolved_boundaries = + resolveBoundaries(boundaries, inner_boundaries); + this->boundaries = resolved_boundaries.first; + this->inner_boundaries = resolved_boundaries.second; + + this->alpha_x = + parseAlphas(alpha_x, this->transport_names, this->n_cols, this->n_rows); + + this->alpha_y = + parseAlphas(alpha_y, this->transport_names, this->n_cols, this->n_rows); +} + +InitialList::DiffusionInit::BoundaryMap +RcppListToBoundaryMap(const std::vector &trans_names, + const Rcpp::List &bound_list, uint32_t n_cols, + uint32_t n_rows) { + InitialList::DiffusionInit::BoundaryMap map; + + for (const auto &name : trans_names) { + const Rcpp::List &conc_list = bound_list[name]; + + tug::Boundary bc(n_rows, n_cols); + + for (const auto &[tug_index, r_init_name] : tug_side_mapping) { + const Rcpp::List &side_list = conc_list[tug_index]; + + const Rcpp::NumericVector &type = side_list["type"]; + const Rcpp::NumericVector &value = side_list["value"]; + + if (type.size() != value.size()) { + throw std::runtime_error( + "Boundary type and value are not the same length"); + } + + for (std::size_t i = 0; i < type.size(); i++) { + if (static_cast(type[i]) == tug::BC_TYPE_CONSTANT) { + bc.setBoundaryElementConstant(static_cast(tug_index), i, + value[i]); + } + } + } + + map[name] = bc.serialize(); + } + + return map; +} +static InitialList::DiffusionInit::InnerBoundaryMap +RcppListToInnerBoundaryMap(const std::vector &trans_names, + const Rcpp::List &inner_bound_list, + std::uint32_t n_cols, std::uint32_t n_rows) { + InitialList::DiffusionInit::InnerBoundaryMap map; + + if (inner_bound_list.size() == 0) { + return map; + } + + for (const auto &name : trans_names) { + const Rcpp::List &conc_list = inner_bound_list[name]; + + std::map, TugType> inner_bc; + + const Rcpp::NumericVector &row = conc_list["row"]; + const Rcpp::NumericVector &col = conc_list["col"]; + const Rcpp::NumericVector &value = conc_list["value"]; + + if (row.size() != col.size() || row.size() != value.size()) { + throw std::runtime_error( + "Inner boundary vectors are not the same length"); + } + + for (std::size_t i = 0; i < row.size(); i++) { + inner_bc[std::make_pair(row[i], col[i])] = value[i]; + } + + map[name] = inner_bc; + } + + return map; +} + +InitialList::DiffusionInit InitialList::getDiffusionInit() const { + DiffusionInit diff_init; + + diff_init.dim = this->dim; + + diff_init.n_cols = this->n_cols; + diff_init.n_rows = this->n_rows; + + diff_init.s_cols = this->s_cols; + diff_init.s_rows = this->s_rows; + + diff_init.constant_cells = this->constant_cells; + diff_init.transport_names = this->transport_names; + + diff_init.boundaries = RcppListToBoundaryMap( + this->transport_names, this->boundaries, this->n_cols, this->n_rows); + diff_init.inner_boundaries = + RcppListToInnerBoundaryMap(this->transport_names, this->inner_boundaries, + this->n_cols, this->n_rows); + diff_init.alpha_x = Field(this->alpha_x); + diff_init.alpha_y = Field(this->alpha_y); + + return diff_init; +} + +} // namespace poet \ No newline at end of file diff --git a/src/Init/GridInit.cpp b/src/Init/GridInit.cpp new file mode 100644 index 000000000..39305dbb2 --- /dev/null +++ b/src/Init/GridInit.cpp @@ -0,0 +1,190 @@ +#include "InitialList.hpp" + +#include "PhreeqcInit.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace poet { + +static Rcpp::NumericMatrix +pqcScriptToGrid(std::unique_ptr &phreeqc, RInside &R) { + PhreeqcInit::PhreeqcMat phreeqc_mat = phreeqc->getPhreeqcMat(); + + // add "id" to the front of the column names + + const std::size_t ncols = phreeqc_mat.names.size(); + const std::size_t nrows = phreeqc_mat.values.size(); + + phreeqc_mat.names.insert(phreeqc_mat.names.begin(), "ID"); + + Rcpp::NumericMatrix mat(nrows, ncols + 1); + + for (std::size_t i = 0; i < nrows; i++) { + mat(i, 0) = phreeqc_mat.ids[i]; + for (std::size_t j = 0; j < ncols; ++j) { + mat(i, j + 1) = phreeqc_mat.values[i][j]; + } + } + + Rcpp::colnames(mat) = Rcpp::wrap(phreeqc_mat.names); + + return mat; +} + +static inline Rcpp::List matToGrid(RInside &R, const Rcpp::NumericMatrix &mat, + const Rcpp::NumericMatrix &grid) { + Rcpp::Function pqc_to_grid_R("pqc_to_grid"); + + return pqc_to_grid_R(mat, grid); +} + +static inline std::map +replaceRawKeywordIDs(std::map raws) { + for (auto &raw : raws) { + std::string &s = raw.second; + // find at beginning of line '*_RAW' followed by a number and change this + // number to 1 + std::regex re(R"((RAW\s+)(\d+))"); + s = std::regex_replace(s, re, "RAW 1"); + } + + return raws; +} + +static inline uint32_t getSolutionCount(std::unique_ptr &phreeqc, + const Rcpp::List &initial_grid) { + PhreeqcInit::ModulesArray mod_array; + Rcpp::Function unique_R("unique"); + + std::vector row_ids = + Rcpp::as>(unique_R(initial_grid["ID"])); + + // std::vector sizes_vec(sizes.begin(), sizes.end()); + + // Rcpp::Function modify_sizes("modify_module_sizes"); + // sizes_vec = Rcpp::as>( + // modify_sizes(sizes_vec, phreeqc_mat, initial_grid)); + + // std::copy(sizes_vec.begin(), sizes_vec.end(), sizes.begin()); + + return phreeqc->getModuleSizes(row_ids)[POET_SOL]; +} + +static std::string readFile(const std::string &path) { + std::string string_rpath(PATH_MAX, '\0'); + + if (realpath(path.c_str(), string_rpath.data()) == nullptr) { + throw std::runtime_error("Failed to resolve the realpath to file " + path); + } + + std::ifstream file(string_rpath); + + if (!file.is_open()) { + throw std::runtime_error("Failed to open file: " + path); + } + + std::stringstream buffer; + buffer << file.rdbuf(); + + return buffer.str(); +} + +void InitialList::prepareGrid(const Rcpp::List &grid_input) { + // parse input values + Rcpp::Function unique_R("unique"); + + std::string script; + std::string database; + + if (grid_input.containsElementNamed( + GRID_MEMBER_STR(GridMembers::PQC_SCRIPT_FILE))) { + + script = readFile(Rcpp::as( + grid_input[GRID_MEMBER_STR(GridMembers::PQC_SCRIPT_FILE)])); + } else { + script = Rcpp::as( + grid_input[GRID_MEMBER_STR(GridMembers::PQC_SCRIPT_STRING)]); + } + + if (grid_input.containsElementNamed( + GRID_MEMBER_STR(GridMembers::PQC_DB_FILE))) { + + database = readFile(Rcpp::as( + grid_input[GRID_MEMBER_STR(GridMembers::PQC_DB_FILE)])); + } else { + database = Rcpp::as( + grid_input[GRID_MEMBER_STR(GridMembers::PQC_DB_STRING)]); + } + + this->database = database; + + Rcpp::NumericMatrix grid_def = + grid_input[GRID_MEMBER_STR(GridMembers::GRID_DEF)]; + Rcpp::NumericVector grid_size = + grid_input[GRID_MEMBER_STR(GridMembers::GRID_SIZE)]; + // Rcpp::NumericVector constant_cells = grid["constant_cells"].; + + this->n_rows = grid_def.nrow(); + this->n_cols = grid_def.ncol(); + + this->s_cols = grid_size[0]; + this->s_rows = grid_size[1]; + + this->dim = n_cols == 1 ? 1 : 2; + + if (this->n_cols > 1 && this->n_rows == 1) { + throw std::runtime_error( + "Dimensions of grid definition does not match the expected format " + "(n_rows > 1, n_cols >= 1)."); + } + + if (this->dim != grid_size.size()) { + throw std::runtime_error( + "Dimensions of grid does not match the dimension of grid definition."); + } + + if (this->s_rows <= 0 || this->s_cols <= 0) { + throw std::runtime_error("Grid size must be positive."); + } + + this->phreeqc = std::make_unique(database, script); + + this->phreeqc_mat = pqcScriptToGrid(phreeqc, R); + this->initial_grid = matToGrid(R, this->phreeqc_mat, grid_def); + + const uint32_t solution_count = getSolutionCount(phreeqc, this->initial_grid); + + std::vector colnames = + Rcpp::as>(this->initial_grid.names()); + + this->transport_names = std::vector( + colnames.begin() + 1, + colnames.begin() + 1 + solution_count); // skip ID + + std::map pqc_raw_dumps; + + pqc_raw_dumps = replaceRawKeywordIDs(phreeqc->raw_dumps()); + + this->pqc_ids = + Rcpp::as>(unique_R(this->initial_grid["ID"])); + + for (const auto &id : this->pqc_ids) { + this->pqc_scripts.push_back(pqc_raw_dumps[id]); + this->pqc_exchanger.push_back(phreeqc->getExchanger(id)); + this->pqc_kinetics.push_back(phreeqc->getKineticsNames(id)); + this->pqc_equilibrium.push_back(phreeqc->getEquilibriumNames(id)); + this->pqc_surface_comps.push_back(phreeqc->getSurfaceCompNames(id)); + this->pqc_surface_charges.push_back(phreeqc->getSurfaceChargeNames(id)); + } +} + +} // namespace poet \ No newline at end of file diff --git a/src/Init/InitialList.cpp b/src/Init/InitialList.cpp new file mode 100644 index 000000000..e4bae7093 --- /dev/null +++ b/src/Init/InitialList.cpp @@ -0,0 +1,135 @@ +#include "InitialList.hpp" +#include "DataStructures/NamedVector.hpp" +#include +#include +#include +#include +#include +#include +#include + +namespace poet { +void InitialList::initializeFromList(const Rcpp::List &setup) { + prepareGrid(setup[grid_key]); + initDiffusion(setup[diffusion_key]); + initChemistry(setup[chemistry_key]); +} + +void InitialList::importList(const Rcpp::List &setup, bool minimal) { + this->dim = Rcpp::as(setup[static_cast(ExportList::GRID_DIM)]); + + Rcpp::NumericVector grid_specs = + setup[static_cast(ExportList::GRID_SPECS)]; + this->n_rows = static_cast(grid_specs[0]); + this->n_cols = static_cast(grid_specs[1]); + + Rcpp::NumericVector spatial = + setup[static_cast(ExportList::GRID_SPATIAL)]; + this->s_rows = spatial[0]; + this->s_cols = spatial[1]; + + this->constant_cells = Rcpp::as>( + setup[static_cast(ExportList::GRID_CONSTANT)]); + + this->porosity = Rcpp::as>( + setup[static_cast(ExportList::GRID_POROSITY)]); + + if (!minimal) { + this->initial_grid = + Rcpp::List(setup[static_cast(ExportList::GRID_INITIAL)]); + } + + this->transport_names = Rcpp::as>( + setup[static_cast(ExportList::DIFFU_TRANSPORT)]); + if (!minimal) { + this->boundaries = + Rcpp::List(setup[static_cast(ExportList::DIFFU_BOUNDARIES)]); + this->inner_boundaries = + Rcpp::List(setup[static_cast(ExportList::DIFFU_INNER_BOUNDARIES)]); + this->alpha_x = + Rcpp::List(setup[static_cast(ExportList::DIFFU_ALPHA_X)]); + this->alpha_y = + Rcpp::List(setup[static_cast(ExportList::DIFFU_ALPHA_Y)]); + } + this->database = + Rcpp::as(setup[static_cast(ExportList::CHEM_DATABASE)]); + this->field_header = Rcpp::as>( + setup[static_cast(ExportList::CHEM_FIELD_HEADER)]); + this->pqc_scripts = Rcpp::as>( + setup[static_cast(ExportList::CHEM_PQC_SCRIPTS)]); + this->pqc_ids = Rcpp::as>( + setup[static_cast(ExportList::CHEM_PQC_IDS)]); + this->pqc_solutions = Rcpp::as>( + setup[static_cast(ExportList::CHEM_PQC_SOLUTIONS)]); + this->pqc_solution_primaries = Rcpp::as>( + setup[static_cast(ExportList::CHEM_PQC_SOLUTION_PRIMARY)]); + this->pqc_exchanger = + Rcpp::List(setup[static_cast(ExportList::CHEM_PQC_EXCHANGER)]); + this->pqc_kinetics = + Rcpp::List(setup[static_cast(ExportList::CHEM_PQC_KINETICS)]); + this->pqc_equilibrium = + Rcpp::List(setup[static_cast(ExportList::CHEM_PQC_EQUILIBRIUM)]); + this->pqc_surface_comps = + Rcpp::List(setup[static_cast(ExportList::CHEM_PQC_SURFACE_COMPS)]); + this->pqc_surface_charges = + Rcpp::List(setup[static_cast(ExportList::CHEM_PQC_SURFACE_CHARGES)]); + + this->dht_species = NamedVector( + setup[static_cast(ExportList::CHEM_DHT_SPECIES)]); + + this->interp_species = Rcpp::as>( + setup[static_cast(ExportList::CHEM_INTERP_SPECIES)]); + + this->chem_hooks = + Rcpp::as(setup[static_cast(ExportList::CHEM_HOOKS)]); +} + +Rcpp::List InitialList::exportList() { + Rcpp::List out(static_cast(ExportList::ENUM_SIZE)); + + out[static_cast(ExportList::GRID_DIM)] = this->dim; + out[static_cast(ExportList::GRID_SPECS)] = + Rcpp::NumericVector::create(this->n_rows, this->n_cols); + out[static_cast(ExportList::GRID_SPATIAL)] = + Rcpp::NumericVector::create(this->s_rows, this->s_cols); + out[static_cast(ExportList::GRID_CONSTANT)] = + Rcpp::wrap(this->constant_cells); + out[static_cast(ExportList::GRID_POROSITY)] = Rcpp::wrap(this->porosity); + out[static_cast(ExportList::GRID_INITIAL)] = this->initial_grid; + + out[static_cast(ExportList::DIFFU_TRANSPORT)] = + Rcpp::wrap(this->transport_names); + out[static_cast(ExportList::DIFFU_BOUNDARIES)] = this->boundaries; + out[static_cast(ExportList::DIFFU_INNER_BOUNDARIES)] = + this->inner_boundaries; + out[static_cast(ExportList::DIFFU_ALPHA_X)] = this->alpha_x; + out[static_cast(ExportList::DIFFU_ALPHA_Y)] = this->alpha_y; + + out[static_cast(ExportList::CHEM_DATABASE)] = Rcpp::wrap(this->database); + out[static_cast(ExportList::CHEM_FIELD_HEADER)] = + Rcpp::wrap(this->field_header); + out[static_cast(ExportList::CHEM_PQC_SCRIPTS)] = + Rcpp::wrap(this->pqc_scripts); + out[static_cast(ExportList::CHEM_PQC_IDS)] = Rcpp::wrap(this->pqc_ids); + out[static_cast(ExportList::CHEM_PQC_SOLUTIONS)] = + Rcpp::wrap(this->pqc_solutions); + out[static_cast(ExportList::CHEM_PQC_SOLUTION_PRIMARY)] = + Rcpp::wrap(this->pqc_solution_primaries); + out[static_cast(ExportList::CHEM_PQC_EXCHANGER)] = + Rcpp::wrap(this->pqc_exchanger); + out[static_cast(ExportList::CHEM_PQC_KINETICS)] = + Rcpp::wrap(this->pqc_kinetics); + out[static_cast(ExportList::CHEM_PQC_EQUILIBRIUM)] = + Rcpp::wrap(this->pqc_equilibrium); + out[static_cast(ExportList::CHEM_PQC_SURFACE_COMPS)] = + Rcpp::wrap(this->pqc_surface_comps); + out[static_cast(ExportList::CHEM_PQC_SURFACE_CHARGES)] = + Rcpp::wrap(this->pqc_surface_charges); + out[static_cast(ExportList::CHEM_DHT_SPECIES)] = this->dht_species; + out[static_cast(ExportList::CHEM_INTERP_SPECIES)] = + Rcpp::wrap(this->interp_species); + out[static_cast(ExportList::CHEM_HOOKS)] = this->chem_hooks; + + return out; +} +} // namespace poet \ No newline at end of file diff --git a/src/Init/InitialList.hpp b/src/Init/InitialList.hpp new file mode 100644 index 000000000..bc9555329 --- /dev/null +++ b/src/Init/InitialList.hpp @@ -0,0 +1,240 @@ +#pragma once + +#include "Base/RInsidePOET.hpp" +#include "DataStructures/NamedVector.hpp" +#include "POETInit.hpp" +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include + +namespace poet { + +using TugType = double; + +class InitialList { +public: + InitialList(RInside &R) : R(R){}; + + void initializeFromList(const Rcpp::List &setup); + + void importList(const Rcpp::List &setup, bool minimal = false); + Rcpp::List exportList(); + + Field getInitialGrid() const { return Field(this->initial_grid); } + +private: + RInside &R; + + enum class ExportList { + GRID_DIM, + GRID_SPECS, + GRID_SPATIAL, + GRID_CONSTANT, + GRID_POROSITY, + GRID_INITIAL, + DIFFU_TRANSPORT, + DIFFU_BOUNDARIES, + DIFFU_INNER_BOUNDARIES, + DIFFU_ALPHA_X, + DIFFU_ALPHA_Y, + CHEM_DATABASE, + CHEM_FIELD_HEADER, + CHEM_PQC_SCRIPTS, + CHEM_PQC_IDS, + CHEM_PQC_SOLUTIONS, + CHEM_PQC_SOLUTION_PRIMARY, + CHEM_PQC_EXCHANGER, + CHEM_PQC_KINETICS, + CHEM_PQC_EQUILIBRIUM, + CHEM_PQC_SURFACE_COMPS, + CHEM_PQC_SURFACE_CHARGES, + CHEM_DHT_SPECIES, + CHEM_INTERP_SPECIES, + CHEM_HOOKS, + ENUM_SIZE + }; + + // Grid members + static constexpr const char *grid_key = "Grid"; + + enum class GridMembers { + PQC_SCRIPT_STRING, + PQC_SCRIPT_FILE, + PQC_DB_STRING, + PQC_DB_FILE, + GRID_DEF, + GRID_SIZE, + CONSTANT_CELLS, + POROSITY, + ENUM_SIZE + }; + + static constexpr std::size_t size_GridMembers = + static_cast(InitialList::GridMembers::ENUM_SIZE); + + static constexpr std::array + GridMembersString = {"pqc_in_string", "pqc_in_file", "pqc_db_string", + "pqc_db_file", "grid_def", "grid_size", + "constant_cells", "porosity"}; + + constexpr const char *GRID_MEMBER_STR(GridMembers member) const { + return GridMembersString[static_cast(member)]; + } + + std::unique_ptr phreeqc; + + void prepareGrid(const Rcpp::List &grid_input); + + std::uint8_t dim{0}; + + std::uint32_t n_cols{0}; + std::uint32_t n_rows{0}; + + double s_cols{0}; + double s_rows{0}; + + std::vector constant_cells; + std::vector porosity; + + Rcpp::List initial_grid; + + // No export + Rcpp::NumericMatrix phreeqc_mat; + +public: + struct DiffusionInit { + using BoundaryMap = + std::map>>; + + using InnerBoundaryMap = + std::map, TugType>>; + + uint8_t dim; + std::uint32_t n_cols; + std::uint32_t n_rows; + + double s_cols; + double s_rows; + + std::vector constant_cells; + + std::vector transport_names; + + BoundaryMap boundaries; + InnerBoundaryMap inner_boundaries; + + Field alpha_x; + Field alpha_y; + }; + + DiffusionInit getDiffusionInit() const; + +private: + // Diffusion members + static constexpr const char *diffusion_key = "Diffusion"; + + enum class DiffusionMembers { + BOUNDARIES, + INNER_BOUNDARIES, + ALPHA_X, + ALPHA_Y, + ENUM_SIZE + }; + + static constexpr std::size_t size_DiffusionMembers = + static_cast(InitialList::DiffusionMembers::ENUM_SIZE); + + static constexpr std::array + DiffusionMembersString = {"boundaries", "inner_boundaries", "alpha_x", + "alpha_y"}; + + constexpr const char *DIFFU_MEMBER_STR(DiffusionMembers member) const { + return DiffusionMembersString[static_cast(member)]; + } + + void initDiffusion(const Rcpp::List &diffusion_input); + std::pair + resolveBoundaries(const Rcpp::List &boundaries_list, + const Rcpp::List &inner_boundaries); + + Rcpp::List boundaries; + Rcpp::List inner_boundaries; + + Rcpp::List alpha_x; + Rcpp::List alpha_y; + + std::vector transport_names; + + // Chemistry Members + static constexpr const char *chemistry_key = "Chemistry"; + + void initChemistry(const Rcpp::List &chem_input); + + std::vector field_header; + + std::string database; + std::vector pqc_scripts; + std::vector pqc_ids; + + std::vector pqc_solutions; + std::vector pqc_solution_primaries; + Rcpp::List pqc_exchanger; + Rcpp::List pqc_kinetics; + Rcpp::List pqc_equilibrium; + Rcpp::List pqc_surface_comps; + Rcpp::List pqc_surface_charges; + + NamedVector dht_species; + + NamedVector interp_species; + + Rcpp::List chem_hooks; + + const std::set hook_name_list{"dht_fill", "dht_fuzz", + "interp_pre", "interp_post"}; + +public: + struct ChemistryHookFunctions { + RHookFunction dht_fill; + RHookFunction> dht_fuzz; + RHookFunction> interp_pre; + RHookFunction interp_post; + }; + + struct ChemistryInit { + uint32_t total_grid_cells; + + // std::vector field_header; + + std::string database; + // std::vector pqc_scripts; + // std::vector pqc_ids; + + std::map pqc_config; + + // std::vector pqc_sol_order; + + NamedVector dht_species; + NamedVector interp_species; + ChemistryHookFunctions hooks; + }; + + ChemistryInit getChemistryInit() const; +}; +} // namespace poet \ No newline at end of file diff --git a/src/Transport/DiffusionModule.cpp b/src/Transport/DiffusionModule.cpp index abbecce95..2ea80564a 100644 --- a/src/Transport/DiffusionModule.cpp +++ b/src/Transport/DiffusionModule.cpp @@ -4,6 +4,9 @@ ** ** Copyright (C) 2018-2022 Marco De Lucia, Max Luebke (GFZ Potsdam) ** +** Copyright (C) 2023-2024 Marco De Lucia (GFZ Potsdam), Max Luebke (University +** of Potsdam) +** ** POET is free software; you can redistribute it and/or modify it under the ** terms of the GNU General Public License as published by the Free Software ** Foundation; either version 2 of the License, or (at your option) any later @@ -19,174 +22,103 @@ */ #include "DiffusionModule.hpp" -#include "../Base/Macros.hpp" -#include -#include +#include "Base/Macros.hpp" +#include "Init/InitialList.hpp" -#include -#include -#include +#include +#include +#include +#include + +#include "tug/Boundary.hpp" +#include "tug/Grid.hpp" +#include "tug/Simulation.hpp" -#include -#include -#include #include #include using namespace poet; -static constexpr double ZERO_MULTIPLICATOR = 10E-14; +static inline std::vector +MatrixToVec(const Eigen::MatrixX &mat) { + std::vector vec(mat.rows() * mat.cols()); -constexpr std::array borders = { - tug::bc::BC_SIDE_LEFT, tug::bc::BC_SIDE_RIGHT, tug::bc::BC_SIDE_TOP, - tug::bc::BC_SIDE_BOTTOM}; - -inline const char *convert_bc_to_config_file(uint8_t in) { - switch (in) { - case tug::bc::BC_SIDE_TOP: - return "N"; - case tug::bc::BC_SIDE_RIGHT: - return "E"; - case tug::bc::BC_SIDE_BOTTOM: - return "S"; - case tug::bc::BC_SIDE_LEFT: - return "W"; - } - return ""; -} - -DiffusionModule::DiffusionModule(const poet::DiffusionParams &diffu_args, - const poet::GridParams &grid_params) - : n_cells_per_prop(grid_params.total_n) { - this->diff_input.setGridCellN(grid_params.n_cells[0], grid_params.n_cells[1]); - this->diff_input.setDomainSize(grid_params.s_cells[0], - grid_params.s_cells[1]); - - this->dim = grid_params.dim; - - this->initialize(diffu_args, grid_params.total_n); -} - -void DiffusionModule::initialize(const poet::DiffusionParams &args, - std::uint32_t n_grid_cells) { - // const poet::DiffusionParams args(this->R); - - // name of props - // - this->prop_names = Rcpp::as>(args.initial_t.names()); - this->prop_count = this->prop_names.size(); - - // initialize field and alphas - this->alpha.reserve(this->prop_count); - - std::vector> initial_values; - - for (uint32_t i = 0; i < this->prop_count; i++) { - // get alphas - we cannot assume alphas are provided in same order as - // initial input - this->alpha.push_back(args.alpha[this->prop_names[i]]); - - const double val = args.initial_t[prop_names[i]]; - std::vector init_val(n_grid_cells, val); - initial_values.push_back(std::move(init_val)); - - if (this->dim == this->DIM_2D) { - tug::bc::BoundaryCondition bc(diff_input.grid.grid_cells[0], - diff_input.grid.grid_cells[1]); - this->bc_vec.push_back(bc); - } else { - tug::bc::BoundaryCondition bc(diff_input.grid.grid_cells[0]); - this->bc_vec.push_back(bc); + for (std::uint32_t i = 0; i < mat.cols(); i++) { + for (std::uint32_t j = 0; j < mat.rows(); j++) { + vec[j * mat.cols() + i] = mat(j, i); } } - this->t_field = Field(n_grid_cells, initial_values, prop_names); - - // apply boundary conditions to each ghost node - uint8_t bc_size = (this->dim == this->DIM_1D ? 2 : 4); - for (uint8_t i = 0; i < bc_size; i++) { - const auto &side = borders[i]; - std::vector vecinj_i = Rcpp::as>( - args.vecinj_index[convert_bc_to_config_file(side)]); - for (int i = 0; i < this->prop_count; i++) { - std::vector bc_vec = args.vecinj[this->prop_names[i]]; - tug::bc::BoundaryCondition &curr_bc = *(this->bc_vec.begin() + i); - for (int j = 0; j < vecinj_i.size(); j++) { - if (vecinj_i[j] == 0) { - continue; - } - if (this->dim == this->DIM_2D) { - curr_bc(side, j) = {tug::bc::BC_TYPE_CONSTANT, - bc_vec[vecinj_i[j] - 1]}; - } - if (this->dim == this->DIM_1D) { - curr_bc(side) = {tug::bc::BC_TYPE_CONSTANT, bc_vec[vecinj_i[j] - 1]}; - } - } - } - } - - // apply inner grid constant cells - // NOTE: opening a scope here for distinguish variable names - if (args.vecinj_inner.size() != 0) { - // apply inner constant cells for every concentration - for (int i = 0; i < this->prop_count; i++) { - std::vector bc_vec = args.vecinj[this->prop_names[i]]; - tug::bc::BoundaryCondition &curr_bc = *(this->bc_vec.begin() + i); - for (int j = 0; j < args.vecinj_inner.size(); j++) { - std::vector inner_tuple = - Rcpp::as>(args.vecinj_inner[j]); - tug::bc::boundary_condition bc = {tug::bc::BC_TYPE_CONSTANT, - bc_vec[inner_tuple[0] - 1]}; - - this->index_constant_cells.push_back(inner_tuple[1]); - uint32_t x = inner_tuple[1]; - uint32_t y = (this->dim == this->DIM_1D ? 0 : inner_tuple[2]); - - curr_bc.setInnerBC(bc, x, y); - } - } - } + return vec; } -void DiffusionModule::simulate(double dt) { - double sim_a_transport, sim_b_transport; +static inline Eigen::MatrixX +VecToMatrix(const std::vector &vec, std::uint32_t n_rows, + std::uint32_t n_cols) { + Eigen::MatrixX mat(n_rows, n_cols); - sim_b_transport = MPI_Wtime(); - - MSG_NOENDL("DiffusionModule::simulate(): Starting diffusion ..."); - std::cout << std::flush; - - std::vector> field_2d = t_field.As2DVector(); - - this->diff_input.setTimestep(dt); - - for (int i = 0; i < field_2d.size(); i++) { - std::vector in_alpha(this->n_cells_per_prop, this->alpha[i]); - this->diff_input.setBoundaryCondition(this->bc_vec[i]); - - if (this->dim == this->DIM_1D) { - tug::diffusion::BTCS_1D(this->diff_input, field_2d[i].data(), - in_alpha.data()); - } else { - tug::diffusion::ADI_2D(this->diff_input, field_2d[i].data(), - in_alpha.data()); + for (std::uint32_t i = 0; i < n_cols; i++) { + for (std::uint32_t j = 0; j < n_rows; j++) { + mat(j, i) = vec[j * n_cols + i]; } } - t_field = field_2d; - - sim_a_transport = MPI_Wtime(); - - transport_t += sim_a_transport - sim_b_transport; - std::cout << " done in " << sim_a_transport - sim_b_transport << "sec" - << std::endl; + return mat; } -void DiffusionModule::end() { - // R["simtime_transport"] = transport_t; - // R.parseEvalQ("profiling$simtime_transport <- simtime_transport"); +// static constexpr double ZERO_MULTIPLICATOR = 10E-14; + +void DiffusionModule::simulate(double requested_dt) { + MSG("Starting diffusion ..."); + const auto start_diffusion_t = std::chrono::high_resolution_clock::now(); + + const auto &n_rows = this->param_list.n_rows; + const auto &n_cols = this->param_list.n_cols; + + tug::Grid grid(param_list.n_rows, param_list.n_cols); + tug::Boundary boundary(grid); + + grid.setDomain(param_list.s_rows, param_list.s_cols); + + tug::Simulation sim(grid, boundary); + sim.setIterations(1); + + for (const auto &sol_name : this->param_list.transport_names) { + auto &species_conc = this->transport_field[sol_name]; + + Eigen::MatrixX conc = VecToMatrix(species_conc, n_rows, n_cols); + Eigen::MatrixX alpha_x = + VecToMatrix(this->param_list.alpha_x[sol_name], n_rows, n_cols); + Eigen::MatrixX alpha_y = + VecToMatrix(this->param_list.alpha_y[sol_name], n_rows, n_cols); + + boundary.deserialize(this->param_list.boundaries[sol_name]); + + if (!this->param_list.inner_boundaries[sol_name].empty()) { + boundary.setInnerBoundaries(this->param_list.inner_boundaries[sol_name]); + } + + grid.setAlpha(alpha_x, alpha_y); + grid.setConcentrations(conc); + + sim.setTimestep(requested_dt); + + sim.run(); + + species_conc = MatrixToVec(grid.getConcentrations()); + } + + const auto end_diffusion_t = std::chrono::high_resolution_clock::now(); + + const auto consumed_time_seconds = + std::chrono::duration_cast>( + end_diffusion_t - start_diffusion_t); + + MSG("Diffusion done in " + std::to_string(consumed_time_seconds.count()) + + " [s]"); + + transport_t += consumed_time_seconds.count(); } double DiffusionModule::getTransportTime() { return this->transport_t; } diff --git a/src/Transport/DiffusionModule.hpp b/src/Transport/DiffusionModule.hpp index a31bc0738..e8a576499 100644 --- a/src/Transport/DiffusionModule.hpp +++ b/src/Transport/DiffusionModule.hpp @@ -21,20 +21,13 @@ #ifndef DIFFUSION_MODULE_H #define DIFFUSION_MODULE_H -#include "../Base/Grid.hpp" -#include "../Base/SimParams.hpp" -#include "../DataStructures/DataStructures.hpp" +#include "DataStructures/Field.hpp" +#include "Init/InitialList.hpp" -#include -#include - -#include -#include -#include -#include -#include +#include namespace poet { + /** * @brief Class describing transport simulation * @@ -53,8 +46,8 @@ public: * * @param R RRuntime object */ - DiffusionModule(const poet::DiffusionParams &diffu_args, - const poet::GridParams &grid_params); + DiffusionModule(const InitialList::DiffusionInit &init_list, Field field) + : param_list(init_list), transport_field(field){}; /** * @brief Run simulation for one iteration @@ -64,14 +57,6 @@ public: */ void simulate(double dt); - /** - * @brief End simulation - * - * All measured timings are distributed to the R runtime - * - */ - void end(); - /** * @brief Get the transport time * @@ -84,39 +69,23 @@ public: * * \return Reference to the diffusion field. */ - Field &getField() { return this->t_field; } + Field &getField() { return this->transport_field; } private: /** * @brief Instance of RRuntime * */ - // RRuntime &R; - enum { DIM_1D = 1, DIM_2D }; + InitialList::DiffusionInit param_list; - void initialize(const poet::DiffusionParams &args, - std::uint32_t n_grid_cells); - - uint8_t dim; - - uint32_t prop_count; - - tug::diffusion::TugInput diff_input; - std::vector alpha; - std::vector index_constant_cells; - std::vector prop_names; - - std::vector bc_vec; - Field t_field; - - uint32_t n_cells_per_prop; + Field transport_field; /** * @brief time spent for transport * */ - double transport_t = 0.f; + double transport_t = 0.; }; } // namespace poet diff --git a/src/initializer.cpp b/src/initializer.cpp new file mode 100644 index 000000000..d9a8afadf --- /dev/null +++ b/src/initializer.cpp @@ -0,0 +1,83 @@ +#include "Init/InitialList.hpp" +#include "poet.hpp" + +#include "Base/argh.hpp" + +#include + +#include +#include +#include +#include + +int main(int argc, char **argv) { + + argh::parser cmdl({"-o", "--output"}); + + cmdl.parse(argc, argv); + + if (cmdl[{"-h", "--help"}] || cmdl.pos_args().size() != 2) { + std::cout << "Usage: " << argv[0] + << " [-o, --output output_file]" + " [-s, --setwd] " + " " + << std::endl; + return EXIT_SUCCESS; + } + + RInside R(argc, argv); + + R.parseEvalQ(init_r_library); + + std::string input_script = cmdl.pos_args()[1]; + std::string normalized_path_script; + std::string in_file_name; + + std::string curr_path = + Rcpp::as(Rcpp::Function("normalizePath")(Rcpp::wrap("."))); + + try { + normalized_path_script = + Rcpp::as(Rcpp::Function("normalizePath")(input_script)); + + in_file_name = Rcpp::as( + Rcpp::Function("basename")(normalized_path_script)); + + std::cout << "Using script " << input_script << std::endl; + Rcpp::Function("source")(input_script); + } catch (Rcpp::exception &e) { + std::cerr << "Error while sourcing " << input_script << e.what() + << std::endl; + return EXIT_FAILURE; + } + + std::string output_file; + + cmdl({"-o", "--output"}, + curr_path + "/" + + in_file_name.substr(0, in_file_name.find_last_of('.')) + ".rds") >> + output_file; + + if (cmdl[{"-s", "--setwd"}]) { + const std::string dir_path = Rcpp::as( + Rcpp::Function("dirname")(normalized_path_script)); + + Rcpp::Function("setwd")(Rcpp::wrap(dir_path)); + } + + Rcpp::List setup = R["setup"]; + + poet::InitialList init(R); + + init.initializeFromList(setup); + + // replace file extension by .rds + Rcpp::Function save("saveRDS"); + + save(init.exportList(), Rcpp::wrap(output_file)); + + std::cout << "Saved result to " << output_file << std::endl; + + // parseGrid(R, grid, results); + return EXIT_SUCCESS; +} \ No newline at end of file diff --git a/src/poet.cpp b/src/poet.cpp index 10c706a0d..48142cbf8 100644 --- a/src/poet.cpp +++ b/src/poet.cpp @@ -4,6 +4,8 @@ ** ** Copyright (C) 2018-2022 Marco De Lucia, Max Luebke (GFZ Potsdam) ** +** Copyright (C) 2023-2024 Max Luebke (University of Potsdam) +** ** POET is free software; you can redistribute it and/or modify it under the ** terms of the GNU General Public License as published by the Free Software ** Foundation; either version 2 of the License, or (at your option) any later @@ -18,123 +20,237 @@ ** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. */ -#include "Base/Grid.hpp" #include "Base/Macros.hpp" #include "Base/RInsidePOET.hpp" -#include "Base/SimParams.hpp" #include "Chemistry/ChemistryModule.hpp" +#include "DataStructures/Field.hpp" +#include "Init/InitialList.hpp" #include "Transport/DiffusionModule.hpp" -#include - +#include #include -#include +#include +#include +#include #include - -#include -#include -#include -#include - +#include #include +#include +#include + +#include "Base/argh.hpp" + +#include using namespace std; using namespace poet; using namespace Rcpp; -poet::ChemistryModule::SingleCMap DFToHashMap(const Rcpp::DataFrame &df) { - std::unordered_map out_map; - vector col_names = Rcpp::as>(df.names()); +static int MY_RANK = 0; - for (const auto &name : col_names) { - double val = df[name.c_str()]; - out_map.insert({name, val}); - } +static std::unique_ptr global_rt_setup; - return out_map; +// we need some layz evaluation, as we can't define the functions before the R +// runtime is initialized +static std::optional master_init_R; +static std::optional master_iteration_end_R; +static std::optional store_setup_R; + +static void init_global_functions(RInside &R) { + R.parseEval(kin_r_library); + master_init_R = Rcpp::Function("master_init"); + master_iteration_end_R = Rcpp::Function("master_iteration_end"); + store_setup_R = Rcpp::Function("StoreSetup"); } // HACK: this is a step back as the order and also the count of fields is // predefined, but it will change in the future -void writeFieldsToR(RInside &R, const Field &trans, const Field &chem) { +// static inline void writeFieldsToR(RInside &R, const Field &trans, +// const Field &chem) { + +// Rcpp::DataFrame t_field(trans.asSEXP()); +// R["TMP"] = t_field; +// R.parseEval("mysetup$state_T <- TMP"); + +// R["TMP"] = chem.asSEXP(); +// R.parseEval("mysetup$state_C <- TMP"); +// } + +enum ParseRet { PARSER_OK, PARSER_ERROR, PARSER_HELP }; + +ParseRet parseInitValues(char **argv, RuntimeParameters ¶ms) { + // initialize argh object + argh::parser cmdl(argv); + + // if user asked for help + if (cmdl[{"help", "h"}]) { + if (MY_RANK == 0) { + MSG("Todo"); + MSG("See README.md for further information."); + } + + return ParseRet::PARSER_HELP; + } + // if positional arguments are missing + if (!cmdl(3)) { + if (MY_RANK == 0) { + ERRMSG("POET needs 3 positional arguments: "); + ERRMSG("1) the R script defining your simulation runtime."); + ERRMSG("2) the initial .rds file generated by poet_init."); + ERRMSG("3) the directory prefix where to save results and profiling"); + } + return ParseRet::PARSER_ERROR; + } + + // parse flags and check for unknown + for (const auto &option : cmdl.flags()) { + if (!(flaglist.find(option) != flaglist.end())) { + if (MY_RANK == 0) { + ERRMSG("Unrecognized option: " + option); + ERRMSG("Make sure to use available options. Exiting!"); + } + return ParseRet::PARSER_ERROR; + } + } + + // parse parameters and check for unknown + for (const auto &option : cmdl.params()) { + if (!(paramlist.find(option.first) != paramlist.end())) { + if (MY_RANK == 0) { + ERRMSG("Unrecognized option: " + option.first); + ERRMSG("Make sure to use available options. Exiting!"); + } + return ParseRet::PARSER_ERROR; + } + } + + params.print_progressbar = cmdl[{"P", "progress"}]; + + /*Parse work package size*/ + cmdl("work-package-size", CHEM_DEFAULT_WORK_PACKAGE_SIZE) >> + params.work_package_size; + + /* Parse DHT arguments */ + params.use_dht = cmdl["dht"]; + params.use_interp = cmdl["interp"]; + // cout << "CPP: DHT is " << ( dht_enabled ? "ON" : "OFF" ) << '\n'; + + cmdl("dht-size", CHEM_DHT_SIZE_PER_PROCESS_MB) >> params.dht_size; + // cout << "CPP: DHT size per process (Byte) = " << dht_size_per_process << + // endl; + + cmdl("dht-snaps", 0) >> params.dht_snaps; + + params.use_interp = cmdl["interp"]; + cmdl("interp-size", 100) >> params.interp_size; + cmdl("interp-min", 5) >> params.interp_min_entries; + cmdl("interp-bucket-entries", 20) >> params.interp_bucket_entries; + + if (MY_RANK == 0) { + // MSG("Complete results storage is " + BOOL_PRINT(simparams.store_result)); + MSG("Work Package Size: " + std::to_string(params.work_package_size)); + MSG("DHT is " + BOOL_PRINT(params.use_dht)); + + if (params.use_dht) { + // MSG("DHT strategy is " + std::to_string(simparams.dht_strategy)); + // MDL: these should be outdated (?) + // MSG("DHT key default digits (ignored if 'signif_vector' is " + // "defined) = " + // << simparams.dht_significant_digits); + // MSG("DHT logarithm before rounding: " + // << (simparams.dht_log ? "ON" : "OFF")); + MSG("DHT size per process (Megabyte) = " + + std::to_string(params.dht_size)); + MSG("DHT save snapshots is " + BOOL_PRINT(params.dht_snaps)); + // MSG("DHT load file is " + chem_params.dht_file); + } + + if (params.use_interp) { + MSG("PHT interpolation enabled: " + BOOL_PRINT(params.use_interp)); + MSG("PHT interp-size = " + std::to_string(params.interp_size)); + MSG("PHT interp-min = " + std::to_string(params.interp_min_entries)); + MSG("PHT interp-bucket-entries = " + + std::to_string(params.interp_bucket_entries)); + } + } + + std::string init_file; + std::string runtime_file; + + cmdl(1) >> runtime_file; + cmdl(2) >> init_file; + cmdl(3) >> params.out_dir; + + // chem_params.dht_outdir = out_dir; + + /* distribute information to R runtime */ + // if local_rank == 0 then master else worker + // R["local_rank"] = MY_RANK; + // assign a char* (string) to 'filesim' + // R["filesim"] = wrap(runtime_file); + // assign a char* (string) to 'fileout' + // R["fileout"] = wrap(out_dir); + // pass the boolean "store_result" to the R process + // R["store_result"] = simparams.store_result; + // // worker count + // R["n_procs"] = simparams.world_size - 1; + // // work package size + // R["work_package_size"] = simparams.wp_size; + // // dht enabled? + // R["dht_enabled"] = chem_params.use_dht; + // // log before rounding? + // R["dht_log"] = simparams.dht_log; + + try { + Rcpp::Function source("source"); + Rcpp::Function readRDS("readRDS"); + + Rcpp::List init_params_ = readRDS(init_file); + params.init_params = init_params_; + + global_rt_setup = std::make_unique(); + *global_rt_setup = source(runtime_file, Rcpp::Named("local", true)); + *global_rt_setup = global_rt_setup->operator[]("value"); + + params.timesteps = + Rcpp::as>(global_rt_setup->operator[]("timesteps")); + + } catch (const std::exception &e) { + ERRMSG("Error while parsing R scripts: " + std::string(e.what())); + return ParseRet::PARSER_ERROR; + } + + return ParseRet::PARSER_OK; +} + +// HACK: this is a step back as the order and also the count of fields is +// predefined, but it will change in the future +void call_master_iter_end(RInside &R, const Field &trans, const Field &chem) { R["TMP"] = Rcpp::wrap(trans.AsVector()); R["TMP_PROPS"] = Rcpp::wrap(trans.GetProps()); - R.parseEval(std::string( - "mysetup$state_T <- setNames(data.frame(matrix(TMP, nrow=" + - std::to_string(trans.GetRequestedVecSize()) + ")), TMP_PROPS)")); + R.parseEval(std::string("state_T <- setNames(data.frame(matrix(TMP, nrow=" + + std::to_string(trans.GetRequestedVecSize()) + + ")), TMP_PROPS)")); R["TMP"] = Rcpp::wrap(chem.AsVector()); R["TMP_PROPS"] = Rcpp::wrap(chem.GetProps()); - R.parseEval(std::string( - "mysetup$state_C <- setNames(data.frame(matrix(TMP, nrow=" + - std::to_string(chem.GetRequestedVecSize()) + ")), TMP_PROPS)")); + R.parseEval(std::string("state_C <- setNames(data.frame(matrix(TMP, nrow=" + + std::to_string(chem.GetRequestedVecSize()) + + ")), TMP_PROPS)")); + R["setup"] = *global_rt_setup; + R.parseEval("setup <- master_iteration_end(setup, state_T, state_C)"); + *global_rt_setup = R["setup"]; } -void set_chem_parameters(poet::ChemistryModule &chem, uint32_t wp_size, - const std::string &database_path) { - chem.SetErrorHandlerMode(1); - chem.SetComponentH2O(false); - chem.SetRebalanceFraction(0.5); - chem.SetRebalanceByCell(true); - chem.UseSolutionDensityVolume(false); - chem.SetPartitionUZSolids(false); +static Rcpp::List RunMasterLoop(RInsidePOET &R, const RuntimeParameters ¶ms, + DiffusionModule &diffusion, + ChemistryModule &chem) { - // Set concentration units - // 1, mg/L; 2, mol/L; 3, kg/kgs - chem.SetUnitsSolution(2); - // 0, mol/L cell; 1, mol/L water; 2 mol/L rock - chem.SetUnitsPPassemblage(1); - // 0, mol/L cell; 1, mol/L water; 2 mol/L rock - chem.SetUnitsExchange(1); - // 0, mol/L cell; 1, mol/L water; 2 mol/L rock - chem.SetUnitsSurface(1); - // 0, mol/L cell; 1, mol/L water; 2 mol/L rock - chem.SetUnitsGasPhase(1); - // 0, mol/L cell; 1, mol/L water; 2 mol/L rock - chem.SetUnitsSSassemblage(1); - // 0, mol/L cell; 1, mol/L water; 2 mol/L rock - chem.SetUnitsKinetics(1); - - // Set representative volume - std::vector rv; - rv.resize(wp_size, 1.0); - chem.SetRepresentativeVolume(rv); - - // Set initial porosity - std::vector por; - por.resize(wp_size, 1); - chem.SetPorosity(por); - - // Set initial saturation - std::vector sat; - sat.resize(wp_size, 1.0); - chem.SetSaturation(sat); - - // Load database - chem.LoadDatabase(database_path); -} - -inline double RunMasterLoop(SimParams ¶ms, RInside &R, - const GridParams &g_params, uint32_t nxyz_master) { - - DiffusionParams d_params{R}; - DiffusionModule diffusion(d_params, g_params); /* Iteration Count is dynamic, retrieving value from R (is only needed by * master for the following loop) */ - uint32_t maxiter = R.parseEval("mysetup$iterations"); + uint32_t maxiter = params.timesteps.size(); - double sim_time = .0; - - ChemistryModule chem(nxyz_master, params.getNumParams().wp_size, maxiter, - params.getChemParams(), MPI_COMM_WORLD); - - set_chem_parameters(chem, nxyz_master, params.getChemParams().database_path); - chem.RunInitFile(params.getChemParams().input_script); - - poet::ChemistryModule::SingleCMap init_df = DFToHashMap(d_params.initial_t); - chem.initializeField(diffusion.getField()); - - if (params.getNumParams().print_progressbar) { + if (params.print_progressbar) { chem.setProgressBarPrintout(true); } @@ -143,47 +259,33 @@ inline double RunMasterLoop(SimParams ¶ms, RInside &R, double dSimTime{0}; for (uint32_t iter = 1; iter < maxiter + 1; iter++) { double start_t = MPI_Wtime(); - uint32_t tick = 0; - // cout << "CPP: Evaluating next time step" << endl; - // R.parseEvalQ("mysetup <- master_iteration_setup(mysetup)"); - double dt = Rcpp::as( - R.parseEval("mysetup$timesteps[" + std::to_string(iter) + "]")); + const double &dt = params.timesteps[iter - 1]; // cout << "CPP: Next time step is " << dt << "[s]" << endl; MSG("Next time step is " + std::to_string(dt) + " [s]"); /* displaying iteration number, with C++ and R iterator */ MSG("Going through iteration " + std::to_string(iter)); - MSG("R's $iter: " + - std::to_string((uint32_t)(R.parseEval("mysetup$iter"))) + - ". Iteration"); /* run transport */ - // TODO: transport to diffusion diffusion.simulate(dt); chem.getField().update(diffusion.getField()); + // chem.getfield().update(diffusion.getfield()); + MSG("Chemistry step"); - chem.SetTimeStep(dt); - chem.RunCells(); - - writeFieldsToR(R, diffusion.getField(), chem.GetField()); - diffusion.getField().update(chem.GetField()); - - R["req_dt"] = dt; - R["simtime"] = (sim_time += dt); - - R.parseEval("mysetup$req_dt <- req_dt"); - R.parseEval("mysetup$simtime <- simtime"); + chem.simulate(dt); // MDL master_iteration_end just writes on disk state_T and // state_C after every iteration if the cmdline option // --ignore-results is not given (and thus the R variable // store_result is TRUE) - R.parseEvalQ("mysetup <- master_iteration_end(setup=mysetup)"); + call_master_iter_end(R, diffusion.getField(), chem.getField()); + + diffusion.getField().update(chem.getField()); MSG("End of *coupling* iteration " + std::to_string(iter) + "/" + std::to_string(maxiter)); @@ -194,175 +296,146 @@ inline double RunMasterLoop(SimParams ¶ms, RInside &R, dSimTime += end_t - start_t; } // END SIMULATION LOOP - R.parseEvalQ("profiling <- list()"); + Rcpp::List chem_profiling; + chem_profiling["simtime"] = chem.GetChemistryTime(); + chem_profiling["loop"] = chem.GetMasterLoopTime(); + chem_profiling["sequential"] = chem.GetMasterSequentialTime(); + chem_profiling["idle_master"] = chem.GetMasterIdleTime(); + chem_profiling["idle_worker"] = Rcpp::wrap(chem.GetWorkerIdleTimings()); + chem_profiling["phreeqc_time"] = Rcpp::wrap(chem.GetWorkerPhreeqcTimings()); - R["simtime_chemistry"] = chem.GetChemistryTime(); - R.parseEvalQ("profiling$simtime_chemistry <- simtime_chemistry"); + Rcpp::List diffusion_profiling; + diffusion_profiling["simtime"] = diffusion.getTransportTime(); - R["chemistry_loop"] = chem.GetMasterLoopTime(); - R.parseEvalQ("profiling$chemistry_loop <- chemistry_loop"); - - R["chemistry_sequential"] = chem.GetMasterSequentialTime(); - R.parseEvalQ("profiling$simtime_sequential <- chemistry_sequential"); - - R["idle_master"] = chem.GetMasterIdleTime(); - R.parseEvalQ("profiling$idle_master <- idle_master"); - - R["idle_worker"] = Rcpp::wrap(chem.GetWorkerIdleTimings()); - R.parseEvalQ("profiling$idle_worker <- idle_worker"); - - R["phreeqc_time"] = Rcpp::wrap(chem.GetWorkerPhreeqcTimings()); - R.parseEvalQ("profiling$phreeqc <- phreeqc_time"); - - R["simtime_transport"] = diffusion.getTransportTime(); - R.parseEvalQ("profiling$simtime_transport <- simtime_transport"); - - // R["phreeqc_count"] = phreeqc_counts; - // R.parseEvalQ("profiling$phreeqc_count <- phreeqc_count"); - - if (params.getChemParams().use_dht) { - R["dht_hits"] = Rcpp::wrap(chem.GetWorkerDHTHits()); - R.parseEvalQ("profiling$dht_hits <- dht_hits"); - R["dht_evictions"] = Rcpp::wrap(chem.GetWorkerDHTEvictions()); - R.parseEvalQ("profiling$dht_evictions <- dht_evictions"); - R["dht_get_time"] = Rcpp::wrap(chem.GetWorkerDHTGetTimings()); - R.parseEvalQ("profiling$dht_get_time <- dht_get_time"); - R["dht_fill_time"] = Rcpp::wrap(chem.GetWorkerDHTFillTimings()); - R.parseEvalQ("profiling$dht_fill_time <- dht_fill_time"); + if (params.use_dht) { + chem_profiling["dht_hits"] = Rcpp::wrap(chem.GetWorkerDHTHits()); + chem_profiling["dht_evictions"] = Rcpp::wrap(chem.GetWorkerDHTEvictions()); + chem_profiling["dht_get_time"] = Rcpp::wrap(chem.GetWorkerDHTGetTimings()); + chem_profiling["dht_fill_time"] = + Rcpp::wrap(chem.GetWorkerDHTFillTimings()); } - if (params.getChemParams().use_interp) { - R["interp_w"] = Rcpp::wrap(chem.GetWorkerInterpolationWriteTimings()); - R.parseEvalQ("profiling$interp_write <- interp_w"); - R["interp_r"] = Rcpp::wrap(chem.GetWorkerInterpolationReadTimings()); - R.parseEvalQ("profiling$interp_read <- interp_r"); - R["interp_g"] = Rcpp::wrap(chem.GetWorkerInterpolationGatherTimings()); - R.parseEvalQ("profiling$interp_gather <- interp_g"); - R["interp_fc"] = + + if (params.use_interp) { + chem_profiling["interp_w"] = + Rcpp::wrap(chem.GetWorkerInterpolationWriteTimings()); + chem_profiling["interp_r"] = + Rcpp::wrap(chem.GetWorkerInterpolationReadTimings()); + chem_profiling["interp_g"] = + Rcpp::wrap(chem.GetWorkerInterpolationGatherTimings()); + chem_profiling["interp_fc"] = Rcpp::wrap(chem.GetWorkerInterpolationFunctionCallTimings()); - R.parseEvalQ("profiling$interp_function_calls <- interp_fc"); - R["interp_calls"] = Rcpp::wrap(chem.GetWorkerInterpolationCalls()); - R.parseEvalQ("profiling$interp_calls <- interp_calls"); - R["interp_cached"] = Rcpp::wrap(chem.GetWorkerPHTCacheHits()); - R.parseEvalQ("profiling$interp_cached <- interp_cached"); + chem_profiling["interp_calls"] = + Rcpp::wrap(chem.GetWorkerInterpolationCalls()); + chem_profiling["interp_cached"] = Rcpp::wrap(chem.GetWorkerPHTCacheHits()); } + Rcpp::List profiling; + profiling["simtime"] = dSimTime; + profiling["chemistry"] = chem_profiling; + profiling["diffusion"] = diffusion_profiling; + chem.MasterLoopBreak(); - diffusion.end(); - return dSimTime; + return profiling; } int main(int argc, char *argv[]) { - double dSimTime, sim_end; - - int world_size, world_rank; + int world_size; MPI_Init(&argc, &argv); - MPI_Comm_size(MPI_COMM_WORLD, &world_size); + { + MPI_Comm_size(MPI_COMM_WORLD, &world_size); + MPI_Comm_rank(MPI_COMM_WORLD, &MY_RANK); - MPI_Comm_rank(MPI_COMM_WORLD, &world_rank); + RInsidePOET &R = RInsidePOET::getInstance(); - RInsidePOET &R = RInsidePOET::getInstance(); - - if (world_rank == 0) { - MSG("Running POET version " + std::string(poet_version)); - } - - if (world_rank > 0) { - { - SimParams params(world_rank, world_size); - int pret = params.parseFromCmdl(argv, R); - - if (pret == poet::PARSER_ERROR) { - MPI_Finalize(); - return EXIT_FAILURE; - } else if (pret == poet::PARSER_HELP) { - MPI_Finalize(); - return EXIT_SUCCESS; - } - - // ChemistryModule worker(nxyz, nxyz, MPI_COMM_WORLD); - ChemistryModule worker = poet::ChemistryModule::createWorker( - MPI_COMM_WORLD, params.getChemParams()); - set_chem_parameters(worker, worker.GetWPSize(), - params.getChemParams().database_path); - - worker.WorkerLoop(); + if (MY_RANK == 0) { + MSG("Running POET version " + std::string(poet_version)); } + RuntimeParameters run_params; + + switch (parseInitValues(argv, run_params)) { + case ParseRet::PARSER_ERROR: + case ParseRet::PARSER_HELP: + MPI_Finalize(); + return 0; + case ParseRet::PARSER_OK: + break; + } + + InitialList init_list(R); + init_list.importList(run_params.init_params, MY_RANK != 0); + + MSG("RInside initialized on process " + std::to_string(MY_RANK)); + + std::cout << std::flush; + MPI_Barrier(MPI_COMM_WORLD); - MSG("finished, cleanup of process " + std::to_string(world_rank)); + ChemistryModule chemistry(run_params.work_package_size, + init_list.getChemistryInit(), MPI_COMM_WORLD); - MPI_Finalize(); + const ChemistryModule::SurrogateSetup surr_setup = { + init_list.getInitialGrid().GetProps(), + run_params.use_dht, + run_params.dht_size, + run_params.use_interp, + run_params.interp_bucket_entries, + run_params.interp_size, + run_params.interp_min_entries}; - return EXIT_SUCCESS; + chemistry.masterEnableSurrogates(surr_setup); + + if (MY_RANK > 0) { + chemistry.WorkerLoop(); + } else { + + init_global_functions(R); + + // R.parseEvalQ("mysetup <- setup"); + // // if (MY_RANK == 0) { // get timestep vector from + // // grid_init function ... // + *global_rt_setup = + master_init_R.value()(*global_rt_setup, run_params.out_dir, + init_list.getInitialGrid().asSEXP()); + + // MDL: store all parameters + // MSG("Calling R Function to store calling parameters"); + // R.parseEvalQ("StoreSetup(setup=mysetup)"); + + MSG("Init done on process with rank " + std::to_string(MY_RANK)); + + // MPI_Barrier(MPI_COMM_WORLD); + + DiffusionModule diffusion(init_list.getDiffusionInit(), + init_list.getInitialGrid()); + + chemistry.masterSetField(init_list.getInitialGrid()); + + Rcpp::List profiling = RunMasterLoop(R, run_params, diffusion, chemistry); + + MSG("finished simulation loop"); + + R["profiling"] = profiling; + R["setup"] = *global_rt_setup; + + string r_vis_code; + r_vis_code = + "saveRDS(profiling, file=paste0(setup$out_dir,'/timings.rds'));"; + R.parseEval(r_vis_code); + + MSG("Done! Results are stored as R objects into <" + run_params.out_dir + + "/timings.rds>"); + } } - /*Loading Dependencies*/ - // TODO: kann raus - R.parseEvalQ(kin_r_library); + MSG("finished, cleanup of process " + std::to_string(MY_RANK)); - SimParams params(world_rank, world_size); - int pret = params.parseFromCmdl(argv, R); - - if (pret == poet::PARSER_ERROR) { - MPI_Finalize(); - return EXIT_FAILURE; - } else if (pret == poet::PARSER_HELP) { - MPI_Finalize(); - return EXIT_SUCCESS; - } - - MSG("RInside initialized on process " + std::to_string(world_rank)); - - R.parseEvalQ("mysetup <- setup"); - // if (world_rank == 0) { // get timestep vector from - // grid_init function ... // - std::string master_init_code = "mysetup <- master_init(setup=setup)"; - R.parseEval(master_init_code); - - GridParams g_params(R); - - params.initVectorParams(R); - - // MDL: store all parameters - if (world_rank == 0) { - MSG("Calling R Function to store calling parameters"); - R.parseEvalQ("StoreSetup(setup=mysetup)"); - } - - if (world_rank == 0) { - MSG("Init done on process with rank " + std::to_string(world_rank)); - } - - // MPI_Barrier(MPI_COMM_WORLD); - - uint32_t nxyz_master = (world_size == 1 ? g_params.total_n : 1); - - dSimTime = RunMasterLoop(params, R, g_params, nxyz_master); - - MSG("finished simulation loop"); - - MSG("start timing profiling"); - - R["simtime"] = dSimTime; - R.parseEvalQ("profiling$simtime <- simtime"); - - string r_vis_code; - r_vis_code = "saveRDS(profiling, file=paste0(fileout,'/timings.rds'));"; - R.parseEval(r_vis_code); - - MSG("Done! Results are stored as R objects into <" + params.getOutDir() + - "/timings.rds>"); - - MPI_Barrier(MPI_COMM_WORLD); - - MSG("finished, cleanup of process " + std::to_string(world_rank)); MPI_Finalize(); - if (world_rank == 0) { + if (MY_RANK == 0) { MSG("done, bye!"); } diff --git a/src/poet.hpp.in b/src/poet.hpp.in index 666fbaa08..f9a86be77 100644 --- a/src/poet.hpp.in +++ b/src/poet.hpp.in @@ -1,10 +1,95 @@ -#ifndef POET_H -#define POET_H +/* +** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of +** Potsdam) +** +** Copyright (C) 2018-2023 Marco De Lucia, Max Luebke (GFZ Potsdam) +** +** Copyright (C) 2023-2024 Max Luebke (University of Potsdam) +** +** POET is free software; you can redistribute it and/or modify it under the +** terms of the GNU General Public License as published by the Free Software +** Foundation; either version 2 of the License, or (at your option) any later +** version. +** +** POET is distributed in the hope that it will be useful, but WITHOUT ANY +** WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +** A PARTICULAR PURPOSE. See the GNU General Public License for more details. +** +** You should have received a copy of the GNU General Public License along with +** this program; if not, write to the Free Software Foundation, Inc., 51 +** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +*/ +#pragma once + +#include +#include #include +#include + +#include + static const char *poet_version = "@POET_VERSION@"; // using the Raw string literal to avoid escaping the quotes -static inline std::string kin_r_library = R"(@R_KIN_LIB@)"; +static const inline std::string kin_r_library = R"(@R_KIN_LIB@)"; -#endif // POET_H +static const inline std::string init_r_library = R"(@R_INIT_LIB@)"; + +static const inline std::string r_runtime_parameters = "mysetup"; + +const std::set flaglist{"ignore-result", "dht", "P", "progress", + "interp"}; +const std::set paramlist{ + "work-package-size", "dht-strategy", "dht-size", "dht-snaps", + "dht-file", "interp-size", "interp-min", "interp-bucket-entries"}; + +constexpr uint32_t CHEM_DEFAULT_WORK_PACKAGE_SIZE = 32; + +constexpr uint32_t CHEM_DHT_SIZE_PER_PROCESS_MB = 1.5E3; + +struct RuntimeParameters { + std::string out_dir; + std::vector timesteps; + + bool print_progressbar; + uint32_t work_package_size; + + Rcpp::List init_params; + + bool use_dht; + std::uint32_t dht_size; + std::uint8_t dht_snaps; + + bool use_interp; + std::uint32_t interp_size; + std::uint32_t interp_min_entries; + std::uint32_t interp_bucket_entries; + + struct ChemistryParams { + // std::string database_path; + // std::string input_script; + + // bool use_dht; + // std::uint64_t dht_size; + // int dht_snaps; + // std::string dht_file; + // std::string dht_outdir; + // NamedVector dht_signifs; + + // bool use_interp; + // std::uint64_t pht_size; + // std::uint32_t pht_max_entries; + // std::uint32_t interp_min_entries; + // NamedVector pht_signifs; + + // struct Chem_Hook_Functions { + // RHookFunction dht_fill; + // RHookFunction> dht_fuzz; + // RHookFunction> interp_pre; + // RHookFunction interp_post; + // } hooks; + + // void initFromR(RInsidePOET &R); + }; +}; diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 884750735..8bc494158 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -2,8 +2,10 @@ file(GLOB test_SRC CONFIGURE_DEPENDS "*.cpp" "*.c") +find_package(doctest REQUIRED) + add_executable(testPOET ${test_SRC}) -target_link_libraries(testPOET doctest poetlib) +target_link_libraries(testPOET doctest::doctest POETLib) target_include_directories(testPOET PRIVATE "${PROJECT_SOURCE_DIR}/src") get_filename_component(TEST_RInsideSourceFile "RInsidePOET_funcs.R" REALPATH) diff --git a/test/RInsidePOET_funcs.R b/test/RInsidePOET_funcs.R index 668b05169..787a42252 100644 --- a/test/RInsidePOET_funcs.R +++ b/test/RInsidePOET_funcs.R @@ -13,10 +13,13 @@ bool_named_vec <- function(input) { } simple_field <- function(field) { + field <- as.data.frame(field, check.names = FALSE) field$Na <- 0 return(field) } extended_field <- function(field, additional) { + field <- as.data.frame(field, check.names = FALSE) + additional <- as.data.frame(additional, check.names = FALSE) return(field + additional) } diff --git a/test/testField.cpp b/test/testField.cpp index acfe8d1a8..0800b4dbd 100644 --- a/test/testField.cpp +++ b/test/testField.cpp @@ -1,15 +1,12 @@ #include -#include -#include #include #include #include -#include #include #include #include -#include +#include #include "testDataStructures.hpp" diff --git a/test/testNamedVector.cpp b/test/testNamedVector.cpp index d5ee908f5..7b86c7496 100644 --- a/test/testNamedVector.cpp +++ b/test/testNamedVector.cpp @@ -1,11 +1,10 @@ #include #include #include -#include #include #include -#include +#include #include "testDataStructures.hpp" diff --git a/util/ci.Dockerfile b/util/ci.Dockerfile deleted file mode 100644 index e795879a5..000000000 --- a/util/ci.Dockerfile +++ /dev/null @@ -1,11 +0,0 @@ -FROM debian - -RUN apt-get update && export DEBIAN_FRONTEND=noninteractive \ - && apt-get install -y \ - cmake \ - git \ - libeigen3-dev \ - libopenmpi-dev \ - r-cran-rcpp \ - r-cran-rinside \ - libssl-dev \ No newline at end of file diff --git a/util/data_evaluation/jlFun_Eval.jl b/util/data_evaluation/jlFun_Eval.jl new file mode 100644 index 000000000..4c95d5ced --- /dev/null +++ b/util/data_evaluation/jlFun_Eval.jl @@ -0,0 +1,31 @@ +using RData, DataFrames, Plots + +# Load all .rds files in given directory +function load_data(directory) + files = readdir(directory) + # data is dictionary with iteration number as key + data = Dict{Int,Any}() + for file in files + if (endswith(file, ".rds") && startswith(file, "iter")) + # extract iteration number (iter_XXX.rds) + iter = parse(Int, split(split(file, "_")[2], ".")[1]) + data[iter] = load(joinpath(directory, file)) + end + end + return data +end + +function spec_to_mat(in_df::DataFrame, spec::Symbol, cols::Int) + specvec = in_df[!, spec] + specmat = transpose(reshape(specvec, cols, :)) + + return specmat +end + +function plot_field(data::AbstractMatrix, log::Bool=false) + if log + heatmap(1:size(data, 2), 1:size(data, 1), log10.(data), c=:viridis) + else + heatmap(1:size(data, 2), 1:size(data, 1), data, c=:viridis) + end +end \ No newline at end of file