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

Compare commits

base: max:main
Branches Tags
max:main
max:naa-naaice
max:naaice-baseline
max:ml/add_real_barite_benchmark
max:control-loop-proto1
max:control-loop-proto2
max:readme
max:wip-selected-output
max:refac/hd5
max:tmp-fgcs
max:serialization-change
max:naaice
max:mdl/fmt
max:ml/field_memory
max:0.4
max:archive/iccs2025
max:archive/hans_ai_surrogate
max:archive/hans_k_means
max:v0.3
max:v0.2.1
max:archive/ns-daos-diffusion-2
max:archive/ns-daos-diffusion-1
max:v0.2
max:v0.1
..
compare: max:archive/ns-daos-diffusion-2
Branches Tags
max:naa-naaice
max:naaice-baseline
max:ml/add_real_barite_benchmark
max:control-loop-proto1
max:control-loop-proto2
max:main
max:readme
max:wip-selected-output
max:refac/hd5
max:tmp-fgcs
max:serialization-change
max:naaice
max:mdl/fmt
max:ml/field_memory
max:0.4
max:archive/iccs2025
max:archive/hans_ai_surrogate
max:archive/hans_k_means
max:v0.3
max:v0.2.1
max:archive/ns-daos-diffusion-2
max:archive/ns-daos-diffusion-1
max:v0.2
max:v0.1

1 Commits
main ... archive/ns

Author SHA1 Message Date
Daos Test
1070fc115b Daos 2023-09-09 01:07:52 +02:00
157 changed files with 17355 additions and 6442 deletions
Show Stats Expand all files Collapse all files

108
.devcontainer/Dockerfile
View File

@ -1,108 +0,0 @@
FROM gcc:11.2.0 AS builder
ENV DEBIAN_FRONTEND=noninteractive
RUN apt-get update \
&& apt-get install -y \
sudo \
git \
ninja-build \
libmpfr-dev \
python3-dev && \
apt-get clean && \
rm -rf /var/lib/apt/lists/*
WORKDIR /tmp
ARG OPENMPI_VERSION=4.1.1
ADD https://download.open-mpi.org/release/open-mpi/v${OPENMPI_VERSION%.*}/openmpi-${OPENMPI_VERSION}.tar.gz /tmp/openmpi.tar.gz
RUN mkdir openmpi && \
tar xf openmpi.tar.gz -C openmpi --strip-components 1 && \
cd openmpi && \
./configure --prefix=/usr/local && \
make -j $(nproc) && \
make install && \
rm -rf /tmp/openmpi tmp/openmpi.tar.gz
ARG CMAKE_VERSION=3.30.5
ADD https://github.com/Kitware/CMake/releases/download/v${CMAKE_VERSION}/cmake-${CMAKE_VERSION}-linux-x86_64.sh /tmp/cmake.sh
RUN bash ./cmake.sh --skip-license --prefix=/usr/local \
&& rm cmake.sh
ARG LAPACK_VERSION=3.12.0
ADD https://github.com/Reference-LAPACK/lapack/archive/refs/tags/v${LAPACK_VERSION}.tar.gz /tmp/lapack.tar.gz
RUN mkdir lapack && \
tar xf lapack.tar.gz -C lapack --strip-components 1 && \
cd lapack && \
mkdir build && \
cd build && \
cmake .. -G Ninja -DBUILD_SHARED_LIBS=ON && \
ninja install && \
rm -rf /tmp/lapack tmp/lapack.tar.gz
ARG R_VERSION=4.4.2
ADD https://cran.r-project.org/src/base/R-${R_VERSION%%.*}/R-${R_VERSION}.tar.gz /tmp/R.tar.gz
RUN mkdir R && \
tar xf R.tar.gz -C R --strip-components 1 && \
cd R && \
./configure --prefix=/usr/local --enable-R-shlib --with-blas --with-lapack && \
make -j $(nproc) && \
make install && \
rm -rf /tmp/R tmp/R.tar.gz
RUN /usr/local/bin/R -q -e "install.packages(c('Rcpp', 'RInside', 'qs'), repos='https://cran.rstudio.com/')"
ARG EIGEN_VERSION=3.4.0
ADD https://gitlab.com/libeigen/eigen/-/archive/${EIGEN_VERSION}/eigen-${EIGEN_VERSION}.tar.bz2 /tmp/eigen.tar.bz2
RUN mkdir eigen && \
tar xf eigen.tar.bz2 -C eigen --strip-components 1 && \
cd eigen && \
mkdir build && \
cd build && \
cmake .. -G Ninja && \
ninja install && \
rm -rf /tmp/eigen tmp/eigen.tar.bz2
ARG GDB_VERSION=15.2
ADD https://ftp.gnu.org/gnu/gdb/gdb-${GDB_VERSION}.tar.xz /tmp/gdb.tar.xz
RUN mkdir gdb && \
tar xf gdb.tar.xz -C gdb --strip-components 1 && \
cd gdb && \
./configure --prefix=/usr/local && \
make -j $(nproc) && \
make install && \
rm -rf /tmp/gdb tmp/gdb.tar.xz
RUN useradd -m -s /bin/bash -G sudo vscode \
&& echo "vscode ALL=(ALL) NOPASSWD:ALL" >> /etc/sudoers
USER vscode
ENV LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH
RUN sudo apt-get update && \
sudo apt-get install -y zsh && \
sudo apt-get clean && \
sudo rm -rf /var/lib/apt/lists/*
RUN sh -c "$(wget -O- https://github.com/deluan/zsh-in-docker/releases/download/v1.2.1/zsh-in-docker.sh)" -- \
-t agnoster \
-p zsh-syntax-highlighting
RUN zsh -c "git clone https://github.com/zsh-users/zsh-syntax-highlighting.git ${ZSH_CUSTOM:-~/.oh-my-zsh/custom}/plugins/zsh-syntax-highlighting"
RUN zsh -c "git clone --depth 1 https://github.com/junegunn/fzf.git ~/.fzf && ~/.fzf/install"
RUN mkdir -p /home/vscode/.config/gdb \
&& echo "set auto-load safe-path /" > /home/vscode/.config/gdb/gdbinit
ENV CMAKE_GENERATOR=Ninja
ENV CMAKE_EXPORT_COMPILE_COMMANDS=ON
WORKDIR /home/vscode

29
.devcontainer/devcontainer.json
View File

@ -1,29 +0,0 @@
// For format details, see https://aka.ms/devcontainer.json. For config options, see the
// README at: https://github.com/devcontainers/templates/tree/main/src/docker-existing-dockerfile
{
"build": {
"dockerfile": "Dockerfile"
},
// Features to add to the dev container. More info: https://containers.dev/features.
// "features": {},
// Use 'forwardPorts' to make a list of ports inside the container available locally.
// "forwardPorts": [],
// Uncomment the next line to run commands after the container is created.
// "postCreateCommand": "cat /etc/os-release",
// Configure tool-specific properties.
"customizations": {
"vscode": {
"extensions": [
"twxs.cmake",
"llvm-vs-code-extensions.vscode-clangd"
]
}
},
// in case you want to push/pull from remote repositories using ssh
"mounts": [
"source=${localEnv:HOME}/.ssh,target=/home/vscode/.ssh,type=bind,consistency=cached",
"source=${localEnv:HOME}/.gitconfig,target=/home/vscode/.gitconfig,type=bind,consistency=cached"
]
// Uncomment to connect as an existing user other than the container default. More info: https://aka.ms/dev-containers-non-root.
// "remoteUser": "devcontainer"
}

3
.gitignore vendored
View File

@ -141,6 +141,3 @@ vignettes/*.pdf
build/
/.cache/
.vscode
.codechecker

151
.gitlab-ci.yml
View File

@ -19,6 +19,7 @@
image: git.gfz-potsdam.de:5000/naaice/poet:ci
stages: # List of stages for jobs, and their order of execution
- build
- release
- test
@ -27,93 +28,89 @@ variables:
SOURCE_ARCHIVE_NAME: 'poet_${CI_COMMIT_TAG}_sources.tar.gz'
CHANGELOG_FILE: 'commit_changelog.md'
test: # This job runs in the build stage, which runs first.
stage: test
build-poet: # This job runs in the build stage, which runs first.
stage: build
script:
- mkdir -p build && cd build
- cmake -DPOET_ENABLE_TESTING=ON -DPOET_PREPROCESS_BENCHS=OFF -DCMAKE_BUILD_TYPE=Release ..
- mkdir build && cd build
- cmake -DPOET_ENABLE_TESTING=ON ..
- make -j$(nproc)
artifacts:
paths:
- build
expire_in: 1 day
test-poet:
stage: test
dependencies:
- build-poet
script:
- cd build
- make -j$(nproc) check
archive-sources: # This job runs in the build stage, which runs first.
image: python:3
stage: release
before_script:
- pip install git-archive-all
- echo ARCHIVE_JOB_ID=${CI_JOB_ID} >> archives.env
script:
- git-archive-all ${SOURCE_ARCHIVE_NAME}
artifacts:
paths:
- ${SOURCE_ARCHIVE_NAME}
expire_in: never
reports:
dotenv: archives.env
rules:
- if: $CI_COMMIT_TAG
release-description:
image: golang:bullseye
stage: release
rules:
- if: $CI_COMMIT_TAG
before_script:
- go install github.com/git-chglog/git-chglog/cmd/git-chglog@v0.15.2
script:
- git-chglog -o ${CHANGELOG_FILE} ${CI_COMMIT_TAG}
artifacts:
paths:
- ${CHANGELOG_FILE}
release-create:
stage: release
image: registry.gitlab.com/gitlab-org/release-cli:latest
rules:
- if: $CI_COMMIT_TAG
script:
- echo "Running release job"
needs:
- job: archive-sources
artifacts: true
- job: release-description
artifacts: true
release:
tag_name: $CI_COMMIT_TAG
name: 'POET $CI_COMMIT_TAG'
description: ${CHANGELOG_FILE}
assets:
links:
- name: '${SOURCE_ARCHIVE_NAME}'
url: 'https://git.gfz-potsdam.de/naaice/poet/-/jobs/${ARCHIVE_JOB_ID}/artifacts/file/${SOURCE_ARCHIVE_NAME}'
pages:
stage: release
before_script:
- apt-get update && apt-get install -y doxygen graphviz
- mkdir {build_pages,public}
- mkdir {build,public}
script:
- pushd build_pages
- pushd build
- cmake .. && make doxygen
- popd && mv build_pages/docs/html/* public/
- popd && mv build/docs/html/* public/
artifacts:
paths:
- public
rules:
- if: $CI_COMMIT_REF_NAME == $CI_DEFAULT_BRANCH || $CI_COMMIT_TAG
push:
stage: release
variables:
GITHUB_REPOSITORY: 'git@github.com:POET-Simulator/POET.git'
before_script:
# I know that there is this file env variable in gitlab, but somehow it does not work for me (still complaining about white spaces ...)
# Therefore, the ssh key is stored as a base64 encoded string
- mkdir -p ~/.ssh && echo $GITHUB_SSH_PRIVATE_KEY | base64 -d > ~/.ssh/id_ed25519 && chmod 0600 ~/.ssh/id_ed25519
- ssh-keyscan github.com >> ~/.ssh/known_hosts
- echo $MIRROR_SCRIPT | base64 -d > mirror.sh && chmod +x mirror.sh
script:
- if [[-d poet.git ]]; then rm -rf poet.git; fi
- git clone --mirror "https://git.gfz-potsdam.de/naaice/poet.git" "poet.git" && cd poet.git
- git push --mirror $GITHUB_REPOSITORY
allow_failure: true
#archive-sources: # This job runs in the build stage, which runs first.
# image: python:3
# stage: release
#
# before_script:
# - pip install git-archive-all
# - echo ARCHIVE_JOB_ID=${CI_JOB_ID} >> archives.env
# script:
# - git-archive-all ${SOURCE_ARCHIVE_NAME}
# artifacts:
# paths:
# - ${SOURCE_ARCHIVE_NAME}
# expire_in: never
# reports:
# dotenv: archives.env
# rules:
# - if: $CI_COMMIT_TAG
#release-description:
# image: golang:bullseye
# stage: release
# rules:
# - if: $CI_COMMIT_TAG
# before_script:
# - go install github.com/git-chglog/git-chglog/cmd/git-chglog@v0.15.2
# script:
# - git-chglog -o ${CHANGELOG_FILE} ${CI_COMMIT_TAG}
# artifacts:
# paths:
# - ${CHANGELOG_FILE}
#
#
#release-create:
# stage: release
# image: registry.gitlab.com/gitlab-org/release-cli:latest
# rules:
# - if: $CI_COMMIT_TAG
# script:
# - echo "Running release job"
# needs:
# - job: archive-sources
# artifacts: true
# - job: release-description
# artifacts: true
# release:
# tag_name: $CI_COMMIT_TAG
# name: 'POET $CI_COMMIT_TAG'
# description: ${CHANGELOG_FILE}
# assets:
# links:
# - name: '${SOURCE_ARCHIVE_NAME}'
# url: 'https://git.gfz-potsdam.de/naaice/poet/-/jobs/${ARCHIVE_JOB_ID}/artifacts/file/${SOURCE_ARCHIVE_NAME}'
- if: $CI_COMMIT_REF_NAME == $CI_DEFAULT_BRANCH

9
.gitmodules vendored
View File

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

51
CITATION.cff
View File

@ -1,51 +0,0 @@
# This CITATION.cff file was generated with cffinit.
# Visit https://bit.ly/cffinit to generate yours today!
cff-version: 1.2.0
title: 'POET: POtsdam rEactive Transport'
message: >-
If you use this software, please cite it using these
metadata.
type: software
authors:
- given-names: Max
family-names: Lübke
email: mluebke@uni-potsdam.de
affiliation: University of Potsdam
orcid: 'https://orcid.org/0009-0008-9773-3038'
- given-names: Marco
family-names: De Lucia
email: delucia@gfz.de
affiliation: GFZ Helmholtz Centre for Geosciences
orcid: 'https://orcid.org/0000-0002-1186-4491'
- given-names: Alexander
family-names: Lindemann
- given-names: Hannes
family-names: Signer
email: signer@uni-potsdam.de
orcid: 'https://orcid.org/0009-0000-3058-8472'
- given-names: Bettina
family-names: Schnor
email: schnor@cs.uni-potsdam.de
affiliation: University of Potsdam
orcid: 'https://orcid.org/0000-0001-7369-8057'
- given-names: Hans
family-names: Straile
identifiers:
- type: doi
value: 10.5194/gmd-14-7391-2021
description: >-
POET (v0.1): speedup of many-core parallel reactive
transport simulations with fast DHT lookups
repository-code: 'https://git.gfz-potsdam.de/naaice/poet'
abstract: >-
Massively parallel reactive transport simulator exploring
acceleration strategies such as embedding of AI/ML and
cache of results in Distributed Hash Tables. Developed in
cooperation with computer scientists of University of
Potsdam.
keywords:
- Reactive Transport
- Geochemistry
- AI/ML Surrogate Modelling
license: GPL-2.0

40
CMake/FindRRuntime.cmake
View File

@ -1,27 +1,30 @@
# prepare R environment (Rcpp + RInside)
find_program(R_EXE "R"
REQUIRED
)
find_program(R_EXE "R")
# search for R executable, R header file and library path
execute_process(
if(R_EXE)
execute_process(
COMMAND ${R_EXE} RHOME
OUTPUT_VARIABLE R_ROOT_DIR
OUTPUT_STRIP_TRAILING_WHITESPACE
)
)
find_path(
find_path(
R_INCLUDE_DIR R.h
HINTS /usr/include /usr/local/include /usr/share ${R_ROOT_DIR}/include
PATH_SUFFIXES R/include R
REQUIRED
)
HINTS ${R_ROOT_DIR}
PATHS /usr/include /usr/local/include /usr/share
PATH_SUFFIXES include/R R/include
)
find_library(
find_library(
R_LIBRARY libR.so
HINTS ${R_ROOT_DIR}/lib
REQUIRED
)
)
else()
message(FATAL_ERROR "No R runtime found!")
endif()
mark_as_advanced(R_INCLUDE_DIR R_LIBRARY R_EXE)
set(R_LIBRARIES ${R_LIBRARY})
set(R_INCLUDE_DIRS ${R_INCLUDE_DIR})
@ -42,6 +45,8 @@ find_path(R_Rcpp_INCLUDE_DIR Rcpp.h
HINTS ${RCPP_PATH}
PATH_SUFFIXES include)
mark_as_advanced(R_Rcpp_INCLUDE_DIR)
list(APPEND R_INCLUDE_DIRS ${R_Rcpp_INCLUDE_DIR})
# find RInside libraries and include path
@ -67,11 +72,16 @@ find_path(R_RInside_INCLUDE_DIR RInside.h
list(APPEND R_LIBRARIES ${R_RInside_LIBRARY})
list(APPEND R_INCLUDE_DIRS ${R_RInside_INCLUDE_DIR})
mark_as_advanced(R_RInside_LIBRARY R_RInside_INCLUDE_DIR)
# putting all together into interface library
add_library(RRuntime INTERFACE IMPORTED)
target_link_libraries(RRuntime INTERFACE ${R_LIBRARIES})
target_include_directories(RRuntime INTERFACE ${R_INCLUDE_DIRS})
set_target_properties(
RRuntime PROPERTIES
INTERFACE_LINK_LIBRARIES "${R_LIBRARIES}"
INTERFACE_INCLUDE_DIRECTORIES "${R_INCLUDE_DIRS}"
)
unset(R_LIBRARIES)
unset(R_INCLUDE_DIRS)

4
CMake/POET_Scripts.cmake
View File

@ -13,7 +13,7 @@ macro(get_POET_version)
WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}
OUTPUT_VARIABLE POET_GIT_VERSION
OUTPUT_STRIP_TRAILING_WHITESPACE)
if(NOT POET_GIT_BRANCH STREQUAL "main")
if(NOT POET_GIT_BRANCH STREQUAL "master")
set(POET_VERSION "${POET_GIT_BRANCH}/${POET_GIT_VERSION}")
else()
set(POET_VERSION "${POET_GIT_VERSION}")
@ -21,7 +21,7 @@ macro(get_POET_version)
elseif(EXISTS ${PROJECT_SOURCE_DIR}/.svn)
file(STRINGS .gitversion POET_VERSION)
else()
set(POET_VERSION "not_versioned")
set(POET_VERSION "0.1")
endif()
message(STATUS "Configuring POET version ${POET_VERSION}")

29
CMakeLists.txt
View File

@ -1,23 +1,14 @@
cmake_minimum_required(VERSION 3.20)
# Version 3.9+ offers new MPI package variables
cmake_minimum_required(VERSION 3.9)
project(POET
LANGUAGES CXX C
DESCRIPTION "A coupled reactive transport simulator")
# specify the C++ standard
set(CMAKE_CXX_STANDARD 20)
set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED True)
set(DEFAULT_BUILD_TYPE "Release")
if(NOT CMAKE_BUILD_TYPE)
message(STATUS "Setting build type to '${DEFAULT_BUILD_TYPE}'.")
set(CMAKE_BUILD_TYPE "${DEFAULT_BUILD_TYPE}" CACHE
STRING "Choose the type of build." FORCE)
# Set the possible values of build type for cmake-gui
set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS
"Debug" "Release" "MinSizeRel" "RelWithDebInfo")
endif()
set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE)
include("CMake/POET_Scripts.cmake")
@ -25,27 +16,27 @@ 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)
add_subdirectory(src)
option(POET_PREPROCESS_BENCHS "Preprocess benchmarks" ON)
if (POET_PREPROCESS_BENCHS)
add_subdirectory(bench)
endif()
add_subdirectory(R_lib)
add_subdirectory(app)
add_subdirectory(bench)
# as tug will also pull in doctest as a dependency
set(TUG_ENABLE_TESTING OFF CACHE BOOL "" FORCE)
add_subdirectory(ext/tug EXCLUDE_FROM_ALL)
add_subdirectory(ext/litephreeqc EXCLUDE_FROM_ALL)
add_subdirectory(ext/phreeqcrm 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()

62
LICENSE
View File

@ -1,8 +1,8 @@
GNU GENERAL PUBLIC LICENSE
Version 2, June 1991
Copyright (C) 1989, 1991 Free Software Foundation, Inc.,
<https://fsf.org/>
Copyright (C) 1989, 1991 Free Software Foundation, Inc., <http://fsf.org/>
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
Everyone is permitted to copy and distribute verbatim copies
of this license document, but changing it is not allowed.
@ -278,61 +278,3 @@ PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
POSSIBILITY OF SUCH DAMAGES.
END OF TERMS AND CONDITIONS
How to Apply These Terms to Your New Programs
If you develop a new program, and you want it to be of the greatest
possible use to the public, the best way to achieve this is to make it
free software which everyone can redistribute and change under these terms.
To do so, attach the following notices to the program. It is safest
to attach them to the start of each source file to most effectively
convey the exclusion of warranty; and each file should have at least
the "copyright" line and a pointer to where the full notice is found.
<one line to give the program's name and a brief idea of what it does.>
Copyright (C) <year> <name of author>
This program 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.
This program 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, see <https://www.gnu.org/licenses/>.
Also add information on how to contact you by electronic and paper mail.
If the program is interactive, make it output a short notice like this
when it starts in an interactive mode:
Gnomovision version 69, Copyright (C) year name of author
Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
This is free software, and you are welcome to redistribute it
under certain conditions; type `show c' for details.
The hypothetical commands `show w' and `show c' should show the appropriate
parts of the General Public License. Of course, the commands you use may
be called something other than `show w' and `show c'; they could even be
mouse-clicks or menu items--whatever suits your program.
You should also get your employer (if you work as a programmer) or your
school, if any, to sign a "copyright disclaimer" for the program, if
necessary. Here is a sample; alter the names:
Yoyodyne, Inc., hereby disclaims all copyright interest in the program
`Gnomovision' (which makes passes at compilers) written by James Hacker.
<signature of Moe Ghoul>, 1 April 1989
Moe Ghoul, President of Vice
This General Public License does not permit incorporating your program into
proprietary programs. If your program is a subroutine library, you may
consider it more useful to permit linking proprietary applications with the
library. If this is what you want to do, use the GNU Lesser General
Public License instead of this License.

310
README.md
View File

@ -1,108 +1,77 @@
<!--
Time-stamp: "Last modified 2023-08-02 13:55:11 mluebke"
-->
# POET
**NOTE: GFZ is migrating its domain from <gfz-potsdam.de> to <gfz.de>.
This should be finalized by the end of 2025. We adopt the NEW domain
in all the links given below. If you encounter 'unreachable address'
try the OLD domain.**
[POET](https://doi.org/10.5281/zenodo.4757913) is a coupled reactive transport
simulator implementing a parallel architecture and a fast, original MPI-based
Distributed Hash Table.
[POET](https://doi.org/10.5281/zenodo.4757913) is a coupled reactive
transport simulator implementing a parallel architecture and a fast,
original MPI-based Distributed Hash Table.
![POET's Coupling Scheme](./docs/POET_scheme.svg)
![POET's Coupling Scheme](./docs/20230720_Scheme_POET_en.svg)
## Parsed code documentiation
A parsed version of POET's documentation can be found at [Gitlab
pages](https://naaice.git-pages.gfz.de/poet).
A parsed version of POET's documentiation can be found at [Gitlab
pages](https://naaice.git-pages.gfz-potsdam.de/poet).
## External Libraries
The following external libraries are shipped with POET:
The following external header library is shipped with POET:
- **CLI11** - <https://github.com/CLIUtils/CLI11>
- **litephreeqc**: IPhreeqc
(<https://github.com/usgs-coupled/iphreeqc>) with patches from
GFZ/UP: <https://git.gfz.de/naaice/litephreeqc>
- **tug** - <https://git.gfz.de/naaice/tug>
- **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
## Installation
### Requirements
To compile POET you need following software to be installed:
To compile POET you need several software to be installed:
- C/C++ compiler (tested with GCC)
- MPI-Implementation (tested with OpenMPI and MVAPICH)
- CMake 3.20+
- Eigen3 3.4+ (required by `tug`)
- *optional*: `doxygen` with `dot` bindings for documentation
- R language and environment including headers or `-dev` packages
(distro dependent)
- R language and environment
- CMake 3.9+
- *optional*: `doxygen` with `dot` bindings for documentiation
The following R packages (and their dependencies) must also be
installed:
The following R libraries must then be installed, which will get the
needed dependencies automatically:
- [Rcpp](https://cran.r-project.org/web/packages/Rcpp/index.html)
- [RInside](https://cran.r-project.org/web/packages/RInside/index.html)
- [qs](https://cran.r-project.org/web/packages/qs/index.html)
- [qs2](https://cran.r-project.org/web/packages/qs2/index.html)
This can be simply achieved by issuing the following commands:
```sh
# start R environment
$ R
# install R dependencies (case sensitive!)
> install.packages(c("Rcpp", "RInside","qs","qs2"))
> q(save="no")
```
### Clone the repository
POET can be anonimously cloned from this repo over https. Make sure to
also download the submodules:
```sh
git clone --recurse-submodules https://git.gfz.de/naaice/poet.git
```
The `--recurse-submodules` option is a shorthand for:
```sh
cd poet
git submodule init && git submodule update
```
### Compiling source code
POET is built with CMake. You can generate Makefiles by running the
usual:
The generation of makefiles is done with CMake. You should be able to generate
Makefiles by running:
```sh
mkdir build && cd build
cmake -DCMAKE_BUILD_TYPE=Release ..
cmake ..
```
This will create the directory `build` and processes the CMake files
and generate Makefiles from it. You're now able to run `make` to start
build process.
If everything went well you'll find the executables at
`build/src/poet`, but it is recommended to install the POET project
If everything went well you'll find the executable at
`build/app/poet`, but it is recommended to install the POET project
structure to a desired `CMAKE_INSTALL_PREFIX` with `make install`.
During the generation of Makefiles, various options can be specified
via `cmake -D <option>=<value> [...]`. Currently, there are the
following available options:
- **POET_DHT_Debug**=_boolean_ - toggles the output of detailed
statistics about DHT usage. Defaults to _OFF_.
- **POET_ENABLE_TESTING**=_boolean_ - enables small set of unit tests
(more to come). Defaults to _OFF_.
- **POET_PHT_ADDITIONAL_INFO**=_boolean_ - enabling the count of
accesses to one PHT bucket. Use with caution, as things will get
slowed down significantly. Defaults to _OFF_.
- **POET_PREPROCESS_BENCHS**=*boolean* - enables the preprocessing of
predefined models/benchmarks. Defaults to *ON*.
- **POET_DHT_Debug**=_boolean_ - toggles the output of detailed statistics about
DHT usage. Defaults to _OFF_.
- **POET_ENABLE_TESTING**=_boolean_ - enables small set of unit tests (more to
come). Defaults to _OFF_.
- **POET_PHT_ADDITIONAL_INFO**=_boolean_ - enabling the count of accesses to one
PHT bucket. Use with caution, as things will get slowed down significantly.
Defaults to _OFF_.
### Example: Build from scratch
@ -115,7 +84,7 @@ follows:
$ R
# install R dependencies
> install.packages(c("Rcpp", "RInside","qs","qs2"))
> install.packages(c("Rcpp", "RInside"))
> q(save="no")
# cd into POET project root
@ -123,7 +92,7 @@ $ cd <POET_dir>
# Build process
$ mkdir build && cd build
$ cmake -DCMAKE_INSTALL_PREFIX=/home/<user>/poet -DCMAKE_BUILD_TYPE=Release ..
$ cmake -DCMAKE_INSTALL_PREFIX=/home/<user>/poet ..
$ make -j<max_numprocs>
$ make install
```
@ -138,57 +107,59 @@ The correspondending directory tree would look like this:
```sh
poet
├── bin
│   ├── poet
│   └── poet_init
│ └── poet
├── R_lib
│ └── kin_r_library.R
└── share
└── poet
└── bench
├── barite
│   ├── barite_200.qs2
│   ├── barite_200_rt.R
│   ├── barite_het.qs2
│   └── barite_het_rt.R
│ ├── barite_interp_eval.R
│ ├── barite.pqi
│ ├── barite.R
│ └── db_barite.dat
├── dolo
│   ├── dolo_inner_large.qs2
│   ├── dolo_inner_large_rt.R
│   ├── dolo_interp.qs2
│   └── dolo_interp_rt.R
│ ├── dolo_diffu_inner_large.R
│ ├── dolo_diffu_inner.R
│ ├── dolo_inner.pqi
│ ├── dolo_interp_long.R
│ └── phreeqc_kin.dat
└── surfex
├── PoetEGU_surfex_500.qs2
└── PoetEGU_surfex_500_rt.R
├── ExBase.pqi
├── ex.R
├── SMILE_2021_11_01_TH.dat
├── SurfExBase.pqi
└── surfex.R
```
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 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.
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.
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.
## Running
Run POET by `mpirun ./poet [OPTIONS] <RUNFILE> <SIMFILE>
<OUTPUT_DIRECTORY>` where:
Run POET by `mpirun ./poet <OPTIONS> <SIMFILE> <OUTPUT_DIRECTORY>`
where:
- **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
- **OPTIONS** - runtime parameters (explained below)
- **SIMFILE** - simulation described as R script (e.g.
`<POET_INSTALL_DIR>/share/poet/bench/dolo/dolo_interp_long.R`)
- **OUTPUT_DIRECTORY** - path, where all output of POET should be stored
### POET command line arguments
### Runtime options
The following parameters can be set:
| Option | Value | Description |
|-----------------------------|--------------|----------------------------------------------------------------------------------|
|-----------------------------|--------------|--------------------------------------------------------------------------------------------------------------------------|
| **--work-package-size=** | _1..n_ | size of work packages (defaults to _5_) |
| **-P, --progress** | | show progress bar |
| **--ai-surrogate** | | activates the AI surrogate chemistry model (defaults to _OFF_) |
| **--ignore-result** | | disables store of simulation resuls |
| **--dht** | | enabling DHT usage (defaults to _OFF_) |
| **--qs** | | store results using qs::qsave() (.qs extension) instead of default qs2 (.qs2) |
| **--rds** | | store results using saveRDS() (.rds extension) instead of default qs2 (.qs2) |
| **--dht-strategy=** | _0-1_ | change DHT strategy. **NOT IMPLEMENTED YET** (Defaults to _0_) |
| **--dht-size=** | _1-n_ | size of DHT per process involved in megabyte (defaults to _1000 MByte_) |
| **--dht-snaps=** | _0-2_ | disable or enable storage of DHT snapshots |
@ -210,106 +181,28 @@ Following values can be set:
### Example: Running from scratch
We will continue the above example and start a simulation with
*barite_het*, which simulation files can be found in
`<INSTALL_DIR>/share/poet/barite/barite_het*`. As transport a
heterogeneous diffusion is used. It's a small 2D grid, 2x5 grid,
simulating 50 time steps with a time step size of 100 seconds. To
start the simulation with 4 processes `cd` into your previously
installed POET-dir `<POET_INSTALL_DIR>/bin` and run:
`dolo_diffu_inner.R`. As transport a simple fixed-coefficient diffusion is used.
It's a 2D, 100x100 grid, simulating 10 time steps. To start the simulation with
4 processes `cd` into your previously installed POET-dir
`<POET_INSTALL_DIR>/bin` and run:
```sh
cp ../share/poet/barite/barite_het* .
mpirun -n 4 ./poet barite_het_rt.R barite_het.qs2 output
mpirun -n 4 ./poet ../share/poet/bench/dolo/dolo_diffu_inner.R/ output
```
After a finished simulation all data generated by POET will be found
in the directory `output`.
You might want to use the DHT to cache previously simulated data and
reuse them in further time-steps. Just append `--dht` to the options
of POET to activate the usage of the DHT. Also, after each iteration a
DHT snapshot shall be produced. This is done by appending the
`--dht-snaps=<value>` option. The resulting call would look like this:
You might want to use the DHT to cache previously simulated data and reuse them
in further time-steps. Just append `--dht` to the options of POET to activate
the usage of the DHT. Also, after each iteration a DHT snapshot shall be
produced. This is done by appending the `--dht-snaps=<value>` option. The
resulting call would look like this:
```sh
mpirun -n 4 ./poet --dht --dht-snaps=2 barite_het_rt.R barite_het.qs2 output
mpirun -n 4 ./poet --dht --dht-snaps=2 ../share/poet/bench/dolo/dolo_diffu_inner.R/ output
```
### Example: Preparing Environment and Running with AI surrogate
To run the AI surrogate, you need to install the R package `keras3`. The
compilation process of POET remains the same as shown above.
In the following code block, the installation process on the Turing Cluster is
shown. `miniconda` is used to create a virtual environment to install
tensorflow/keras. Please adapt the installation process to your needs.
<!-- Start an R interactive session and install the required packages: -->
```sh
# First, install the required R packages
R -e "install.packages('keras3', repos='https://cloud.r-project.org/')"
# manually create a virtual environment to install keras/python using conda,
# as this is somehow broken on the Turing Cluster when using the `keras::install_keras()` function
cd poet
# create a virtual environment in the .ai directory with python 3.11
conda create -p ./.ai python=3.11
conda activate ./.ai
# install tensorflow and keras
pip install keras tensorflow[and-cuda]
# add conda's python path to the R environment
# make sure to have the conda environment activated
echo -e "RETICULATE_PYTHON=$(which python)\n" >> ~/.Renviron
```
After setup the R environment, recompile POET and you're ready to run the AI
surrogate.
```sh
cd <installation_dir>/bin
# copy the benchmark files to the installation directory
cp <project_root_dir>/bench/barite/{barite_50ai*,db_barite.dat,barite.pqi} .
# preprocess the benchmark
./poet_init barite_50ai.R
# run POET with AI surrogate and GPU utilization
srun --gres=gpu -N 1 -n 12 ./poet --ai-surrogate barite_50ai_rt.R barite_50ai.qs2 output
```
Keep in mind that the AI surrogate is currently not stable or might also not
produce any valid predictions.
## 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.de/naaice/poet/-/wikis/Initialization).
We try to keep the document up-to-date. However, if you encounter
missing information or need help, please get in touch with us via the
issue tracker or E-Mail.
`poet_init` can be used as follows:
```sh
./poet_init [-o, --output output_file] [-s, --setwd] <script.R>
```
where:
- **output** - name of the output file (defaults to the input file
name with the extension `.qs2`)
- **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
@ -318,44 +211,3 @@ important information from the OpenMPI Man Page:
For example, on platforms that support it, the clock_gettime()
function will be used to obtain a monotonic clock value with whatever
precision is supported on that platform (e.g., nanoseconds).
## Additional functions for the AI surrogate
The AI surrogate can be activated for any benchmark and is by default
initiated as a sequential keras model with three hidden layer of depth
48, 96, 24 with relu activation and adam optimizer. All functions in
`ai_surrogate_model.R` can be overridden by adding custom definitions
via an R file in the input script. This is done by adding the path to
this file in the input script. Simply add the path as an element
called `ai_surrogate_input_script` to the `chemistry_setup` list.
Please use the global variable `ai_surrogate_base_path` as a base path
when relative filepaths are used in custom funtions.
**There is currently no default implementation to determine the
validity of predicted values.** This means, that every input script
must include an R source file with a custom function
`validate_predictions(predictors, prediction)`. Examples for custom
functions can be found for the barite_200 benchmark
The functions can be defined as follows:
`validate_predictions(predictors, prediction)`: Returns a boolean
index vector that signals for each row in the predictions if the
values are considered valid. Can eg. be implemented as a mass balance
threshold between the predictors and the prediction.
`initiate_model()`: Returns a keras model. Can be used to load
pretrained models.
`preprocess(df, backtransform = FALSE, outputs = FALSE)`: Returns the
scaled/transformed/backtransformed dataframe. The `backtransform` flag
signals if the current processing step is applied to data that's
assumed to be scaled and expects backtransformed values. The `outputs`
flag signals if the current processing step is applied to the output
or tatget of the model. This can be used to eg. skip these processing
steps and only scale the model input.
`training_step (model, predictor, target, validity)`: Trains the model
after each iteration. `validity` is the bool index vector given by
`validate_predictions` and can eg. be used to only train on values
that have not been valid predictions.

1
R_lib/CMakeLists.txt Normal file
View File

@ -0,0 +1 @@
install(FILES kin_r_library.R DESTINATION R_lib)

75
R_lib/ai_surrogate_model.R
View File

@ -1,75 +0,0 @@
## This file contains default function implementations for the ai surrogate.
## To load pretrained models, use pre-/postprocessing or change hyperparameters
## it is recommended to override these functions with custom implementations via
## the input script. The path to the R-file containing the functions mus be set
## in the variable "ai_surrogate_input_script". See the barite_200.R file as an
## example and the general README for more information.
require(keras3)
require(tensorflow)
initiate_model <- function() {
hidden_layers <- c(48, 96, 24)
activation <- "relu"
loss <- "mean_squared_error"
input_length <- length(ai_surrogate_species)
output_length <- length(ai_surrogate_species)
## Creates a new sequential model from scratch
model <- keras_model_sequential()
## Input layer defined by input data shape
model %>% layer_dense(units = input_length,
activation = activation,
input_shape = input_length,
dtype = "float32")
for (layer_size in hidden_layers) {
model %>% layer_dense(units = layer_size,
activation = activation,
dtype = "float32")
}
## Output data defined by output data shape
model %>% layer_dense(units = output_length,
activation = activation,
dtype = "float32")
model %>% compile(loss = loss,
optimizer = "adam")
return(model)
}
gpu_info <- function() {
msgm(tf_gpu_configured())
}
prediction_step <- function(model, predictors) {
prediction <- predict(model, as.matrix(predictors))
colnames(prediction) <- colnames(predictors)
return(as.data.frame(prediction))
}
preprocess <- function(df, backtransform = FALSE, outputs = FALSE) {
return(df)
}
postprocess <- function(df, backtransform = TRUE, outputs = TRUE) {
return(df)
}
set_valid_predictions <- function(temp_field, prediction, validity) {
temp_field[validity == 1, ] <- prediction[validity == 1, ]
return(temp_field)
}
training_step <- function(model, predictor, target, validity) {
msgm("Training:")
x <- as.matrix(predictor)
y <- as.matrix(target[colnames(x)])
model %>% fit(x, y)
model %>% save_model_tf(paste0(out_dir, "/current_model.keras"))
}

112
R_lib/init_r_lib.R
View File

@ -1,112 +0,0 @@
### Copyright (C) 2018-2024 Marco De Lucia, Max Luebke (GFZ Potsdam, 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.
##' @param pqc_mat matrix, containing IDs and PHREEQC outputs
##' @param grid matrix, zonation referring to pqc_mat$ID
##' @return a data.frame
# pqc_to_grid <- function(pqc_mat, grid) {
# # Convert the input DataFrame to a matrix
# pqc_mat <- as.matrix(pqc_mat)
# # Flatten the matrix into a vector
# id_vector <- as.integer(t(grid))
# # Find the matching rows in the matrix
# row_indices <- match(id_vector, pqc_mat[, "ID"])
# # Extract the matching rows from pqc_mat to size of grid matrix
# result_mat <- pqc_mat[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)
# }
##' @param pqc_mat matrix, containing IDs and PHREEQC outputs
##' @param grid matrix, zonation referring to pqc_mat$ID
##' @return a data.frame
pqc_to_grid <- function(pqc_mat, grid) {
# Convert the input DataFrame to a matrix
pqc_mat <- as.matrix(pqc_mat)
# Flatten the matrix into a vector
id_vector <- as.integer(t(grid))
# Find the matching rows in the matrix
row_indices <- match(id_vector, pqc_mat[, "ID"])
# Extract the matching rows from pqc_mat to size of grid matrix
result_mat <- pqc_mat[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)
}
##' @param pqc_mat matrix,
##' @param transport_spec column name of species in pqc_mat
##' @param id
##' @title
##' @return
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)
}
##' @title
##' @param init_grid
##' @param new_names
##' @return
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)
}

362
R_lib/kin_r_library.R
View File

@ -1,3 +1,5 @@
## Time-stamp: "Last modified 2023-08-15 11:58:23 delucia"
### 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
@ -13,57 +15,61 @@
### this program; if not, write to the Free Software Foundation, Inc., 51
### Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
master_init <- function(setup, out_dir, init_field) {
## Simple function to check file extension. It is needed to check if
## the GridFile is SUM (MUFITS format) or rds/RData
FileExt <- function(x) {
pos <- regexpr("\\.([[:alnum:]]+)$", x)
ifelse(pos > -1L, substring(x, pos + 1L), "")
}
master_init <- function(setup) {
msgm("Process with rank 0 reading GRID properties")
## Setup the directory where we will store the results
if (!dir.exists(out_dir)) {
dir.create(out_dir)
msgm("created directory ", out_dir)
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 ", out_dir, " already exists, I will overwrite!")
msgm("dir ", fileout, " already exists, I will overwrite!")
}
if (is.null(setup$store_result)) {
if (!exists("store_result")) {
msgm("store_result doesn't exist!")
} else {
msgm("store_result is ", setup$store_result)
msgm("store_result is ", store_result)
}
} else {
}
setup$iter <- 1
setup$maxiter <- setup$iterations
setup$timesteps <- setup$timesteps
setup$maxiter <- length(setup$timesteps)
setup$iterations <- setup$maxiter
setup$simulation_time <- 0
dgts <- as.integer(ceiling(log10(setup$maxiter)))
## string format to use in sprintf
fmt <- paste0("%0", dgts, "d")
if (is.null(setup[["store_result"]])) {
setup$store_result <- TRUE
}
if (setup$store_result) {
init_field_out <- paste0(out_dir, "/iter_", sprintf(fmt = fmt, 0), ".", setup$out_ext)
init_field <- data.frame(init_field, check.names = FALSE)
SaveRObj(x = init_field, path = 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, state_T, state_C) {
master_iteration_end <- function(setup) {
iter <- setup$iter
# print(iter)
## max digits for iterations
dgts <- as.integer(ceiling(log10(setup$maxiter + 1)))
dgts <- as.integer(ceiling(log10(setup$iterations + 1)))
## string format to use in sprintf
fmt <- paste0("%0", dgts, "d")
@ -71,141 +77,217 @@ master_iteration_end <- function(setup, state_T, state_C) {
## comprised in setup$out_save
if (setup$store_result) {
if (iter %in% setup$out_save) {
nameout <- paste0(setup$out_dir, "/iter_", sprintf(fmt = fmt, iter), ".", setup$out_ext)
state_T <- data.frame(state_T, check.names = FALSE)
state_C <- data.frame(state_C, check.names = FALSE)
ai_surrogate_info <- list(
prediction_time = if (exists("ai_prediction_time")) as.integer(ai_prediction_time) else NULL,
training_time = if (exists("ai_training_time")) as.integer(ai_training_time) else NULL,
valid_predictions = if (exists("validity_vector")) validity_vector else NULL
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)
)
SaveRObj(x = list(
T = state_T,
C = state_C,
simtime = as.integer(setup$simulation_time),
totaltime = as.integer(totaltime),
ai_surrogate_info = ai_surrogate_info
), path = nameout)
saveRDS(list(
T = setup$state_T, C = setup$state_C,
simtime = as.integer(setup$simtime),
tr_info = info
), file = nameout)
msgm("results stored in <", nameout, ">")
}
}
## Add last time step to simulation time
setup$simulation_time <- setup$simulation_time + setup$timesteps[iter]
## msgm("done iteration", iter, "/", length(setup$timesteps))
msgm("done iteration", iter, "/", setup$maxiter)
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(...) {
prefix <- paste0("R: ")
if (local_rank == 0) {
fname <- as.list(sys.call(-1))[[1]]
prefix <- paste0("R: ", fname, " ::")
cat(paste(prefix, ..., "\n"))
}
invisible()
}
## Function called by master R process to store on disk all relevant
## parameters for the simulation
StoreSetup <- function(setup) {
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,
## phase = setup$phase,
## density = setup$density,
## dt_differ = setup$dt_differ,
## prolong = setup$prolong,
## maxiter = setup$maxiter,
## saved_iter = setup$iter_output,
## out_save = setup$out_save )
to_store$Transport <- setup$diffusion
## to_store$Chemistry <- list(
## nprocs = n_procs,
## wp_size = work_package_size,
## base = setup$base,
## first = setup$first,
## init = setup$initsim,
## db = db,
## kin = setup$kin,
## ann = setup$ann)
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
}
to_store$Cmdline <- commandArgs(trailingOnly=FALSE)
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)))
}
## Handler to read R objs from binary files using either builtin
## readRDS(), qs::qread() or qs2::qs_read() based on file extension
ReadRObj <- function(path) {
## code borrowed from tools::file_ext()
pos <- regexpr("\\.([[:alnum:]]+)$", path)
extension <- ifelse(pos > -1L, substring(path, pos + 1L), "")
switch(extension,
rds = readRDS(path),
qs = qs::qread(path),
qs2 = qs2::qs_read(path)
)
}
## Handler to store R objs to binary files using either builtin
## saveRDS() or qs::qsave() based on file extension
SaveRObj <- function(x, path) {
## msgm("Storing to", path)
## code borrowed from tools::file_ext()
pos <- regexpr("\\.([[:alnum:]]+)$", path)
extension <- ifelse(pos > -1L, substring(path, pos + 1L), "")
switch(extension,
rds = saveRDS(object = x, file = path),
qs = qs::qsave(x = x, file = path),
qs2 = qs2::qs_save(object = x, file = path)
)
}
######## Old relic code
## ## Function called by master R process to store on disk all relevant
## ## parameters for the simulation
## 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,
## ## phase = setup$phase,
## ## density = setup$density,
## ## dt_differ = setup$dt_differ,
## ## prolong = setup$prolong,
## ## maxiter = setup$maxiter,
## ## saved_iter = setup$iter_output,
## ## out_save = setup$out_save )
## to_store$Transport <- setup$diffusion
## ## to_store$Chemistry <- list(
## ## nprocs = n_procs,
## ## wp_size = work_package_size,
## ## base = setup$base,
## ## first = setup$first,
## ## init = setup$initsim,
## ## db = db,
## ## kin = setup$kin,
## ## ann = setup$ann)
## 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"))
## }

7
app/CMakeLists.txt Normal file
View File

@ -0,0 +1,7 @@
configure_file(poet.h.in poet.h)
add_executable(poet poet.cpp)
target_include_directories(poet PUBLIC "${CMAKE_CURRENT_BINARY_DIR}")
target_link_libraries(poet PUBLIC poet_lib MPI::MPI_CXX)
install(TARGETS poet DESTINATION bin)

369
app/poet.cpp Normal file
View File

@ -0,0 +1,369 @@
/*
** 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.
*/
#include <Rcpp.h>
#include <cstdint>
#include <cstdlib>
#include <poet/ChemistryModule.hpp>
#include <poet/DiffusionModule.hpp>
#include <poet/Grid.hpp>
#include <poet/Macros.hpp>
#include <poet/RInsidePOET.hpp>
#include <poet/SimParams.hpp>
#include <cstring>
#include <iostream>
#include <string>
#include <vector>
#include <mpi.h>
#include <poet.h>
using namespace std;
using namespace poet;
using namespace Rcpp;
poet::ChemistryModule::SingleCMap DFToHashMap(const Rcpp::DataFrame &df) {
std::unordered_map<std::string, double> out_map;
vector<string> col_names = Rcpp::as<vector<string>>(df.names());
for (const auto &name : col_names) {
double val = df[name.c_str()];
out_map.insert({name, val});
}
return out_map;
}
// 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) {
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["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)"));
}
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);
// Set concentration units
// 1, mg/L; 2, mol/L; 3, kg/kgs
chem.SetUnitsSolution(2);
// 0, mol/L cell; 1, mol/L water; 2 mol/L rock
chem.SetUnitsPPassemblage(1);
// 0, mol/L cell; 1, mol/L water; 2 mol/L rock
chem.SetUnitsExchange(1);
// 0, mol/L cell; 1, mol/L water; 2 mol/L rock
chem.SetUnitsSurface(1);
// 0, mol/L cell; 1, mol/L water; 2 mol/L rock
chem.SetUnitsGasPhase(1);
// 0, mol/L cell; 1, mol/L water; 2 mol/L rock
chem.SetUnitsSSassemblage(1);
// 0, mol/L cell; 1, mol/L water; 2 mol/L rock
chem.SetUnitsKinetics(1);
// Set representative volume
std::vector<double> rv;
rv.resize(wp_size, 1.0);
chem.SetRepresentativeVolume(rv);
// Set initial porosity
std::vector<double> por;
por.resize(wp_size, 1);
chem.SetPorosity(por);
// Set initial saturation
std::vector<double> sat;
sat.resize(wp_size, 1.0);
chem.SetSaturation(sat);
// Load database
chem.LoadDatabase(database_path);
}
inline double RunMasterLoop(SimParams &params, RInside &R,
const GridParams &g_params, uint32_t nxyz_master) {
DiffusionParams d_params{R};
DiffusionModule diffusion(d_params, g_params);
/* Iteration Count is dynamic, retrieving value from R (is only needed by
* master for the following loop) */
uint32_t maxiter = R.parseEval("mysetup$iterations");
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) {
chem.setProgressBarPrintout(true);
}
/* SIMULATION LOOP */
double dSimTime{0};
for (uint32_t iter = 1; iter < maxiter + 1; iter++) {
double start_t = MPI_Wtime();
uint32_t tick = 0;
// cout << "CPP: Evaluating next time step" << endl;
// R.parseEvalQ("mysetup <- master_iteration_setup(mysetup)");
double dt = Rcpp::as<double>(
R.parseEval("mysetup$timesteps[" + std::to_string(iter) + "]"));
// 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());
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");
// 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)");
MSG("End of *coupling* iteration " + std::to_string(iter) + "/" +
std::to_string(maxiter));
MSG();
// MPI_Barrier(MPI_COMM_WORLD);
double end_t = MPI_Wtime();
dSimTime += end_t - start_t;
} // END SIMULATION LOOP
R.parseEvalQ("profiling <- list()");
R["simtime_chemistry"] = chem.GetChemistryTime();
R.parseEvalQ("profiling$simtime_chemistry <- simtime_chemistry");
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.getChemParams().use_interp) {
R["interp_w"] = Rcpp::wrap(chem.GetWorkerInterpolationWriteTimings());
R.parseEvalQ("profiling$interp_write <- interp_w");
R["interp_r"] = Rcpp::wrap(chem.GetWorkerInterpolationReadTimings());
R.parseEvalQ("profiling$interp_read <- interp_r");
R["interp_g"] = Rcpp::wrap(chem.GetWorkerInterpolationGatherTimings());
R.parseEvalQ("profiling$interp_gather <- interp_g");
R["interp_fc"] =
Rcpp::wrap(chem.GetWorkerInterpolationFunctionCallTimings());
R.parseEvalQ("profiling$interp_function_calls <- interp_fc");
R["interp_calls"] = Rcpp::wrap(chem.GetWorkerInterpolationCalls());
R.parseEvalQ("profiling$interp_calls <- interp_calls");
R["interp_cached"] = Rcpp::wrap(chem.GetWorkerPHTCacheHits());
R.parseEvalQ("profiling$interp_cached <- interp_cached");
}
chem.MasterLoopBreak();
diffusion.end();
return dSimTime;
}
int main(int argc, char *argv[]) {
double dSimTime, sim_end;
int world_size, world_rank;
MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD, &world_size);
MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);
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();
}
MPI_Barrier(MPI_COMM_WORLD);
MSG("finished, cleanup of process " + std::to_string(world_rank));
MPI_Finalize();
return EXIT_SUCCESS;
}
/*Loading Dependencies*/
// TODO: kann raus
std::string r_load_dependencies = "source('../R_lib/kin_r_library.R');";
R.parseEvalQ(r_load_dependencies);
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) {
MSG("done, bye!");
}
exit(0);
}

11
app/poet.h.in Normal file
View File

@ -0,0 +1,11 @@
#ifndef POET_H
#define POET_H
#include "poet/ChemistryModule.hpp"
#include <Rcpp.h>
const char *poet_version = "@POET_VERSION@";
const char *CHEMISTRY_MODULE_NAME = "state_C";
#endif // POET_H

4
apps/CMakeLists.txt
View File

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

3
apps/initializer/main.cpp
View File

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

43
bench/CMakeLists.txt
View File

@ -1,43 +1,4 @@
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})
set(OUT_FILE_EXT ${OUT_FILE}.qs2)
add_custom_command(
OUTPUT ${OUT_FILE_EXT}
COMMAND $<TARGET_FILE:poet_init> -n ${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_EXT})
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)
add_subdirectory(DAOSBenchmark)

7
bench/DAOSBenchmark/CMakeLists.txt Normal file
View File

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

205
bench/DAOSBenchmark/dolo_diffu_edge.R Normal file
View File

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

28
bench/DAOSBenchmark/dolo_inner.pqi Normal file
View File

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

1307
bench/DAOSBenchmark/phreeqc_kin.dat Normal file
View File

File diff suppressed because it is too large Load Diff

26
bench/barite/CMakeLists.txt
View File

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

14
bench/barite/README.org
View File

@ -18,19 +18,9 @@ mpirun -np 4 ./poet --interp barite_interp_eval.R barite_results
* List of Files
- =barite_het.R=: POET input script with homogeneous zones for a 5x2 simulation
- =barite.R=: POET input script for a 20x20 simulation grid
- =barite_interp_eval.R=: POET input script for a 400x200 simulation
grid
- =barite_200.R=: POET input script for a 200x200 simulation
grid
- =barite_200ai_surrogate_input_script.R=: Defines the ai surrogate functions
to load a pretrained model and apply min-max-feature scaling on the model inputs
and target. Prediction validity is assessed with a threshold of 3e-5 on the mass
balance of Ba and Sr.
- =barite_200min_max_bounds=: Minimum and maximum values from 50 iterations of the
barite_200 benchmark. Used for feature scaling in the ai surrogate.
- =barite_200model_min_max.keras=: A sequential keras model that has been trained
on 50 iterations of the barite_200 benchmark with min-max-scaled inputs
and targets/outputs.
- =db_barite.dat=: PHREEQC database containing the kinetic expressions
for barite and celestite, stripped down from =phreeqc.dat=
- =barite.pqi=: PHREEQC input script defining the chemical system

147
bench/barite/barite.R Normal file
View File

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

53
bench/barite/barite.pqi
View File

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

59
bench/barite/barite_200.R
View File

@ -1,59 +0,0 @@
cols <- 200
rows <- 200
s_cols <- 1
s_rows <- 1
grid_def <- matrix(2, nrow = rows, ncol = cols)
# 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
)
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
diffusion_setup <- list(
boundaries = list(
"W" = bound_def,
"N" = bound_def
),
alpha_x = homogenous_alpha,
alpha_y = homogenous_alpha
)
dht_species <- c(
"H" = 3,
"O" = 3,
"Charge" = 6,
"Ba" = 6,
"Cl" = 6,
"S" = 6,
"Sr" = 6,
"Barite" = 5,
"Celestite" = 5
)
chemistry_setup <- list(
dht_species = dht_species,
ai_surrogate_input_script = "./barite_200ai_surrogate_input_script.R"
)
# 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
)

7
bench/barite/barite_200_rt.R
View File

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

48
bench/barite/barite_200ai_surrogate_input_script.R
View File

@ -1,48 +0,0 @@
## load a pretrained model from tensorflow file
## Use the global variable "ai_surrogate_base_path" when using file paths
## relative to the input script
initiate_model <- function() {
init_model <- normalizePath(paste0(ai_surrogate_base_path,
"model_min_max_float64.keras"))
return(load_model_tf(init_model))
}
scale_min_max <- function(x, min, max, backtransform) {
if (backtransform) {
return((x * (max - min)) + min)
} else {
return((x - min) / (max - min))
}
}
preprocess <- function(df, backtransform = FALSE, outputs = FALSE) {
minmax_file <- normalizePath(paste0(ai_surrogate_base_path,
"min_max_bounds.rds"))
global_minmax <- readRDS(minmax_file)
for (column in colnames(df)) {
df[column] <- lapply(df[column],
scale_min_max,
global_minmax$min[column],
global_minmax$max[column],
backtransform)
}
return(df)
}
mass_balance <- function(predictors, prediction) {
dBa <- abs(prediction$Ba + prediction$Barite -
predictors$Ba - predictors$Barite)
dSr <- abs(prediction$Sr + prediction$Celestite -
predictors$Sr - predictors$Celestite)
return(dBa + dSr)
}
validate_predictions <- function(predictors, prediction) {
epsilon <- 3e-5
mb <- mass_balance(predictors, prediction)
msgm("Mass balance mean:", mean(mb))
msgm("Mass balance variance:", var(mb))
msgm("Rows where mass balance meets threshold", epsilon, ":",
sum(mb < epsilon))
return(mb < epsilon)
}

60
bench/barite/barite_50ai.R
View File

@ -1,60 +0,0 @@
## Time-stamp: "Last modified 2024-05-30 13:34:14 delucia"
cols <- 50
rows <- 50
s_cols <- 0.25
s_rows <- 0.25
grid_def <- matrix(2, nrow = rows, ncol = cols)
# 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
)
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-8
diffusion_setup <- list(
boundaries = list(
"W" = bound_def,
"N" = bound_def
),
alpha_x = homogenous_alpha,
alpha_y = homogenous_alpha
)
dht_species <- c(
"H" = 3,
"O" = 3,
"Charge" = 3,
"Ba" = 6,
"Cl" = 6,
"S" = 6,
"Sr" = 6,
"Barite" = 5,
"Celestite" = 5
)
chemistry_setup <- list(
dht_species = dht_species,
ai_surrogate_input_script = "./barite_50ai_surr_mdl.R"
)
# 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
)

BIN
bench/barite/barite_50ai_all.keras
View File

Binary file not shown.

9
bench/barite/barite_50ai_rt.R
View File

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

90
bench/barite/barite_50ai_surr_mdl.R
View File

@ -1,90 +0,0 @@
## Time-stamp: "Last modified 2024-05-30 13:27:06 delucia"
## load a pretrained model from tensorflow file
## Use the global variable "ai_surrogate_base_path" when using file paths
## relative to the input script
initiate_model <- function() {
require(keras3)
require(tensorflow)
init_model <- normalizePath(paste0(ai_surrogate_base_path,
"barite_50ai_all.keras"))
Model <- keras3::load_model(init_model)
msgm("Loaded model:")
print(str(Model))
return(Model)
}
scale_min_max <- function(x, min, max, backtransform) {
if (backtransform) {
return((x * (max - min)) + min)
} else {
return((x - min) / (max - min))
}
}
minmax <- list(min = c(H = 111.012433592824, O = 55.5062185549492, Charge = -3.1028354471876e-08,
Ba = 1.87312878574393e-141, Cl = 0, `S(6)` = 4.24227510643685e-07,
Sr = 0.00049382996130541, Barite = 0.000999542409828586, Celestite = 0.244801877115968),
max = c(H = 111.012433679682, O = 55.5087003521685, Charge = 5.27666636082035e-07,
Ba = 0.0908849779513762, Cl = 0.195697626449355, `S(6)` = 0.000620774752665846,
Sr = 0.0558680070692722, Barite = 0.756779139057097, Celestite = 1.00075422160624
))
preprocess <- function(df) {
if (!is.data.frame(df))
df <- as.data.frame(df, check.names = FALSE)
as.data.frame(lapply(colnames(df),
function(x) scale_min_max(x=df[x],
min=minmax$min[x],
max=minmax$max[x],
backtransform=FALSE)),
check.names = FALSE)
}
postprocess <- function(df) {
if (!is.data.frame(df))
df <- as.data.frame(df, check.names = FALSE)
as.data.frame(lapply(colnames(df),
function(x) scale_min_max(x=df[x],
min=minmax$min[x],
max=minmax$max[x],
backtransform=TRUE)),
check.names = FALSE)
}
mass_balance <- function(predictors, prediction) {
dBa <- abs(prediction$Ba + prediction$Barite -
predictors$Ba - predictors$Barite)
dSr <- abs(prediction$Sr + prediction$Celestite -
predictors$Sr - predictors$Celestite)
return(dBa + dSr)
}
validate_predictions <- function(predictors, prediction) {
epsilon <- 1E-7
mb <- mass_balance(predictors, prediction)
msgm("Mass balance mean:", mean(mb))
msgm("Mass balance variance:", var(mb))
ret <- mb < epsilon
msgm("Rows where mass balance meets threshold", epsilon, ":",
sum(ret))
return(ret)
}
training_step <- function(model, predictor, target, validity) {
msgm("Starting incremental training:")
## x <- as.matrix(predictor)
## y <- as.matrix(target[colnames(x)])
history <- model %>% keras3::fit(x = data.matrix(predictor),
y = data.matrix(target),
epochs = 10, verbose=1)
keras3::save_model(model,
filepath = paste0(out_dir, "/current_model.keras"),
overwrite=TRUE)
return(model)
}

32
bench/barite/barite_het.R
View File

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

80
bench/barite/barite_het.pqi
View File

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

4
bench/barite/barite_het_rt.R
View File

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

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

@ -0,0 +1,151 @@
## Time-stamp: "Last modified 2023-07-21 15:04:49 mluebke"
database <- normalizePath("../share/poet/bench/barite/db_barite.dat")
input_script <- normalizePath("../share/poet/bench/barite/barite.pqi")
#################################################################
## Section 1 ##
## Grid initialization ##
#################################################################
n <- 400
m <- 200
types <- c("scratch", "phreeqc", "rds")
init_cell <- list(
"H" = 110.0124,
"O" = 55.5087,
"Charge" = -1.217E-09,
"Ba" = 1.E-10,
"Cl" = 2.E-10,
"S" = 6.205E-4,
"Sr" = 6.205E-4,
"Barite" = 0.001,
"Celestite" = 1
)
grid <- list(
n_cells = c(n, m),
s_cells = c(n / 10, m / 10),
type = types[1],
init_cell = as.data.frame(init_cell, check.names = FALSE),
props = names(init_cell),
database = database,
input_script = input_script
)
##################################################################
## Section 2 ##
## Diffusion parameters and boundary conditions ##
##################################################################
## initial conditions
init_diffu <- list(
# "H" = 110.0124,
"H" = 0.00000028904,
# "O" = 55.5087,
"O" = 0.000000165205,
# "Charge" = -1.217E-09,
"Charge" = -3.337E-08,
"Ba" = 1.E-10,
"Cl" = 1.E-10,
"S(6)" = 6.205E-4,
"Sr" = 6.205E-4
)
injection_diff <- list(
list(
# "H" = 111.0124,
"H" = 0.0000002890408,
# "O" = 55.50622,
"O" = 0.00002014464,
# "Charge" = -3.337E-08,
"Charge" = -3.337000004885E-08,
"Ba" = 0.1,
"Cl" = 0.2,
"S(6)" = 0,
"Sr" = 0
)
)
## diffusion coefficients
alpha_diffu <- c(
"H" = 1E-06,
"O" = 1E-06,
"Charge" = 1E-06,
"Ba" = 1E-06,
"Cl" = 1E-06,
"S(6)" = 1E-06,
"Sr" = 1E-06
)
vecinj_inner <- list(
l1 = c(1, floor(n / 2), floor(m / 2))
## l2 = c(2,80,80),
## l3 = c(2,60,80)
)
boundary <- list(
# "N" = rep(1, n),
"N" = rep(0, n),
"E" = rep(0, n),
"S" = rep(0, n),
"W" = rep(0, n)
)
diffu_list <- names(alpha_diffu)
vecinj <- do.call(rbind.data.frame, injection_diff)
names(vecinj) <- names(init_diffu)
diffusion <- list(
init = as.data.frame(init_diffu, check.names = FALSE),
vecinj = vecinj,
vecinj_inner = vecinj_inner,
vecinj_index = boundary,
alpha = alpha_diffu
)
#################################################################
## Section 3 ##
## Chemistry module (Phreeqc) ##
#################################################################
## # Needed when using DHT
dht_species <- c(
"H" = 10,
"O" = 10,
"Charge" = 3,
"Ba" = 5,
"Cl" = 5,
"S(6)" = 5,
"Sr" = 5
)
chemistry <- list(
database = database,
input_script = input_script,
dht_species = dht_species
)
#################################################################
## Section 4 ##
## Putting all those things together ##
#################################################################
iterations <- 200
dt <- 250
setup <- list(
grid = grid,
diffusion = diffusion,
chemistry = chemistry,
iterations = iterations,
timesteps = rep(dt, iterations),
store_result = TRUE,
out_save = seq(1, iterations)
)

BIN
bench/barite/min_max_bounds.rds
View File

Binary file not shown.

BIN
bench/barite/model_min_max_float64.keras
View File

Binary file not shown.

25
bench/dolo/CMakeLists.txt
View File

@ -1,18 +1,9 @@
set(bench_files
dolo_inner_large.R
dolo_interp.R
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(runtime_files
dolo_inner_large_rt.R
dolo_interp_rt.R
)
ADD_BENCH_TARGET(
dolo_bench
bench_files
runtime_files
"dolo"
)
add_dependencies(${BENCHTARGET} dolo_bench)

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

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

9
bench/dolo/README.org
View File

@ -18,13 +18,16 @@ mpirun -np 4 ./poet --dht --interp dolo_interp_long.R dolo_interp_long_res
* List of Files
- =dolo_interp.R=: POET input script for a 400x200 simulation
- =dolo_diffu_inner.R=: POET input script for a 100x100 simulation
grid
- =dolo_interp_long.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

54
bench/dolo/dol.pqi
View File

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

190
bench/dolo/dolo_diffu_inner.R Normal file
View File

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

190
bench/dolo/dolo_diffu_inner_large.R Normal file
View File

@ -0,0 +1,190 @@
## 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(2, 1),
type = types[1]
)
##################################################################
## Section 2 ##
## Diffusion parameters and boundary conditions ##
##################################################################
## initial conditions
init_diffu <- list(
"H" = 0.000211313883539788,
"O" = 0.00398302904424952,
"Charge" = -5.0822e-19,
"C(4)" = 1.2279E-4,
"Ca" = 1.2279E-4,
"Cl" = 0,
"Mg" = 0
)
## diffusion coefficients
alpha_diffu <- c(
"H" = 1E-6,
"O" = 1E-6,
"Charge" = 1E-6,
"C(4)" = 1E-6,
"Ca" = 1E-6,
"Cl" = 1E-6,
"Mg" = 1E-6
)
## list of boundary conditions/inner nodes
vecinj_diffu <- list(
list(
"H" = 0.0001540445,
"O" = 0.002148006,
"Charge" = 1.90431e-16,
"C(4)" = 0,
"Ca" = 0,
"Cl" = 0.002,
"Mg" = 0.001
),
list(
"H" = 0.0001610193,
"O" = 0.002386934,
"Charge" = 1.90431e-16,
"C(4)" = 0,
"Ca" = 0.0,
"Cl" = 0.004,
"Mg" = 0.002
)
)
vecinj_inner <- list(
l1 = c(1, 400, 200),
l2 = c(2, 1400, 800),
l3 = c(2, 1600, 800)
)
boundary <- list(
# "N" = c(1, rep(0, n-1)),
"N" = rep(0, n),
"E" = rep(0, m),
"S" = rep(0, n),
"W" = rep(0, m)
)
diffu_list <- names(alpha_diffu)
vecinj <- do.call(rbind.data.frame, vecinj_diffu)
names(vecinj) <- names(init_diffu)
diffusion <- list(
init = as.data.frame(init_diffu, check.names = FALSE),
vecinj = vecinj,
vecinj_inner = vecinj_inner,
vecinj_index = boundary,
alpha = alpha_diffu
)
#################################################################
## Section 3 ##
## Chemistry module (Phreeqc) ##
#################################################################
## # Needed when using DHT
dht_species <- c(
"H" = 10,
"O" = 10,
"Charge" = 3,
"C(4)" = 5,
"Ca" = 5,
"Cl" = 5,
"Mg" = 5,
"Calcite" = 5,
"Dolomite" = 5
)
check_sign_cal_dol_dht <- function(old, new) {
if ((old["Calcite"] == 0) != (new["Calcite"] == 0)) {
return(TRUE)
}
if ((old["Dolomite"] == 0) != (new["Dolomite"] == 0)) {
return(TRUE)
}
return(FALSE)
}
fuzz_input_dht_keys <- function(input) {
return(input[names(dht_species)])
}
check_sign_cal_dol_interp <- function(to_interp, data_set) {
data_set <- as.data.frame(do.call(rbind, data_set), check.names = FALSE, optional = TRUE)
names(data_set) <- names(dht_species)
cal <- (data_set$Calcite == 0) == (to_interp["Calcite"] == 0)
dol <- (data_set$Dolomite == 0) == (to_interp["Dolomite"] == 0)
cal_dol_same_sig <- cal == dol
return(rev(which(!cal_dol_same_sig)))
}
check_neg_cal_dol <- function(result) {
neg_sign <- (result["Calcite"] <- 0) || (result["Dolomite"] < 0)
return(any(neg_sign))
}
hooks <- list(
dht_fill = check_sign_cal_dol_dht,
dht_fuzz = fuzz_input_dht_keys,
interp_pre_func = check_sign_cal_dol_interp,
interp_post_func = check_neg_cal_dol
)
chemistry <- list(
database = database,
input_script = input_script,
dht_species = dht_species,
hooks = hooks
)
#################################################################
## Section 4 ##
## Putting all those things together ##
#################################################################
iterations <- 500
dt <- 50
setup <- list(
grid = grid,
diffusion = diffusion,
chemistry = chemistry,
iterations = iterations,
timesteps = rep(dt, iterations),
store_result = TRUE,
out_save = seq(5, iterations, by = 5)
)

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

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

BIN
bench/dolo/dolo_inner.rds
View File

Binary file not shown.

115
bench/dolo/dolo_inner_large.R
View File

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

10
bench/dolo/dolo_inner_large_rt.R
View File

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

131
bench/dolo/dolo_interp.R
View File

@ -1,131 +0,0 @@
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" = 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" = 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" = 3,
"Ca" = 3,
"Mg" = 2,
"Calcite" = 2,
"Dolomite" = 2
)
chemistry_setup <- list(
dht_species = c(
"H" = 3,
"O" = 3,
"Charge" = 3,
"C" = 6,
"Ca" = 6,
"Cl" = 3,
"Mg" = 5,
"Calcite" = 4,
"Dolomite" = 4
),
pht_species = pht_species,
hooks = list(
dht_fill = check_sign_cal_dol_dht,
# dht_fuzz = fuzz_input_dht_keys,
interp_pre = check_sign_cal_dol_interp,
interp_post = check_neg_cal_dol
)
)
# Define a setup list for simulation configuration
setup <- list(
Grid = grid_setup, # Parameters related to the grid structure
Diffusion = diffusion_setup, # Parameters related to the diffusion process
Chemistry = chemistry_setup # Parameters related to the chemistry process
)

204
bench/dolo/dolo_interp_long.R Normal file
View File

@ -0,0 +1,204 @@
## 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(any(neg_sign))
}
hooks <- list(
dht_fill = check_sign_cal_dol_dht,
dht_fuzz = fuzz_input_dht_keys,
interp_pre_func = check_sign_cal_dol_interp,
interp_post_func = check_neg_cal_dol
)
chemistry <- list(
database = database,
input_script = input_script,
dht_species = dht_species,
hooks = hooks
# pht_species = pht_species
)
#################################################################
## Section 4 ##
## Putting all those things together ##
#################################################################
iterations <- 20000
dt <- 200
setup <- list(
grid = grid,
diffusion = diffusion,
chemistry = chemistry,
iterations = iterations,
timesteps = rep(dt, iterations),
store_result = TRUE,
out_save = c(1, seq(50, iterations, by = 50))
)

10
bench/dolo/dolo_interp_rt.R
View File

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

BIN
bench/fgcs/20241211_README.pdf
View File

Binary file not shown.

102
bench/fgcs/20241211_README.tex
View File

@ -1,102 +0,0 @@
% Created 2024-12-11 mer 23:24
% Intended LaTeX compiler: pdflatex
\documentclass[a4paper, 9pt]{article}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{graphicx}
\usepackage{longtable}
\usepackage{wrapfig}
\usepackage{rotating}
\usepackage[normalem]{ulem}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{capt-of}
\usepackage{hyperref}
\usepackage{fullpage}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{charter}
\usepackage{listings}
\lstloadlanguages{R}
\author{MDL <delucia@gfz.de>}
\date{2024-12-11}
\title{A \texttt{barite}-based benchmark for FGCS interpolation paper}
\begin{document}
\maketitle
\section{Description}
\label{sec:org739879a}
\begin{itemize}
\item \texttt{barite\_fgcs\_2.R}: POET input script with circular
"crystals" on a 200x200 nodes grid
\item \(\alpha\): isotropic 10\textsuperscript{-5}
m\textsuperscript{2}/s outside of the crystals,
10\textsuperscript{-7} inside
\item 200 iterations, dt = 1000
\item \texttt{barite\_fgcs\_2.pqi}: PHREEQC input, 4 SOLUTIONS
(basically the same as in \texttt{barite} benchmark):
\begin{enumerate}
\item Equilibrium with Celestite, no mineral \(Rightarrow\)
\item Equilibrium with Celestite, KINETICS Celestite (1 mol) and
Barite (0 mol)
\item Injection of 0.1 BaCl2 from NW corner
\item Injection of 0.2 BaCl2 from SE corner
\end{enumerate}
\item \texttt{db\_barite.dat}: PHREEQC database containing the kinetic
expressions for barite and celestite, stripped down from
\texttt{phreeqc.dat}
\end{itemize}
\begin{figure}[htbp]
\centering
\includegraphics[width=0.48\textwidth]{./fgcs_Celestite_init.pdf}
\includegraphics[width=0.48\textwidth]{./fgcs_Barite_200.pdf}
\caption{\textbf{Left:} Initial distribution of Celestite
"crystals". \textbf{Right:} precipitated Barite}
\end{figure}
\section{Interpolation}
\label{sec:org2a09431}
Using the following parametrization:
\begin{lstlisting}
dht_species <- c("H" = 7,
"O" = 7,
"Ba" = 7,
"Cl" = 7,
"S(6)" = 7,
"Sr" = 7,
"Barite" = 4,
"Celestite" = 4)
pht_species <- c("Ba" = 4,
"Cl" = 3,
"S(6)" = 3,
"Sr" = 3,
"Barite" = 2,
"Celestite" = 2 )
\end{lstlisting}
Runtime goes from 1800 to 600 s (21 CPUs) but there are "suspect"
errors especially in O and H, where "suspect" means some values appear
to be multiplied by 2:
\begin{figure}[htbp]
\centering
\includegraphics[width=0.9\textwidth]{./fgcs_interp_1.pdf}
\caption{Scatterplots reference vs interpolated after 1 coupling
iteration}
\end{figure}
\end{document}
%%% Local Variables:
%%% mode: LaTeX
%%% TeX-master: t
%%% End:

90
bench/fgcs/EvalFGCS.R
View File

@ -1,90 +0,0 @@
## Time-stamp: "Last modified 2024-12-11 23:21:25 delucia"
library(PoetUtils)
library(viridis)
res <- ReadPOETSims("./res_fgcs2_96/")
pp <- PlotField(res$iter_200$C$Barite, rows = 200, cols = 200, contour = FALSE,
nlevels=12, palette=terrain.colors)
cairo_pdf("fgcs_Celestite_init.pdf", family="serif")
par(mar=c(0,0,0,0))
pp <- PlotField((res$iter_000$Celestite), rows = 200, cols = 200,
contour = FALSE, breaks=c(-0.5,0.5,1.5),
palette = grey.colors, plot.axes = FALSE, scale = FALSE,
main="Initial Celestite crystals")
dev.off()
cairo_pdf("fgcs_Ba_init.pdf", family="serif")
par(mar=c(0,0,0,0))
pp <- PlotField(log10(res$iter_001$C$Cl), rows = 200, cols = 200,
contour = FALSE,
palette = terrain.colors, plot.axes = FALSE, scale = FALSE,
main="log10(Ba)")
dev.off()
pp <- PlotField(log10(res$iter_002$C$Ba), rows = 200, cols = 200,
contour = FALSE, palette = viridis, rev.palette = FALSE,
main = "log10(Ba) after 5 iterations")
pp <- PlotField(log10(res$iter_200$C$`S(6)`), rows = 200, cols = 200, contour = FALSE)
str(res$iter_00)
res$iter_178$C$Barite
pp <- res$iter_043$C$Barite
breaks <- pretty(pp, n = 5)
br <- c(0, 0.0005, 0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1)
pp <- PlotField(res$iter_200$C$Barite, rows = 200, cols = 200, contour = FALSE,
breaks = br, palette=terrain.colors)
cairo_pdf("fgcs_Barite_200.pdf", family="serif")
pp <- PlotField(log10(res$iter_200$C$Barite), rows = 200, cols = 200,
contour = FALSE, palette = terrain.colors, plot.axes = FALSE,
rev.palette = FALSE, main = "log10(Barite) after 200 iter")
dev.off()
ref <- ReadPOETSims("./res_fgcs_2_ref")
rei <- ReadPOETSims("./res_fgcs_2_interp1/")
timref <- ReadRObj("./res_fgcs_2_ref/timings.qs")
timint <- ReadRObj("./res_fgcs_2_interp1/timings.qs")
timref
timint
wch <- c("H","O", "Ba", "Sr","Cl", "S(6)")
rf <- data.matrix(ref$iter_001$C[, wch])
r1 <- data.matrix(rei$iter_001$C[, wch])
r1[is.nan(r1)] <- NA
rf[is.nan(rf)] <- NA
cairo_pdf("fgcs_interp_1.pdf", family="serif", width = 10, height = 7)
PlotScatter(rf, r1, which = wch, labs = c("ref", "interp"), cols = 3, log="", las = 1, pch=4)
dev.off()
head(r1)
head(rf)
rf$O
r1$O

2
bench/fgcs/README.org
View File

@ -1,2 +0,0 @@
* Refer to the LaTeX file (and pdf) for more information

105
bench/fgcs/barite_fgcs_2.R
View File

@ -1,105 +0,0 @@
## Time-stamp: "Last modified 2024-12-11 16:08:11 delucia"
cols <- 1000
rows <- 1000
dim_cols <- 50
dim_rows <- 50
ncirc <- 20 ## number of crystals
rmax <- cols / 10 ## max radius (in nodes)
set.seed(22933)
centers <- cbind(sample(seq_len(cols), ncirc), sample(seq_len(rows), ncirc))
radii <- sample(seq_len(rmax), ncirc, replace = TRUE)
mi <- matrix(rep(seq_len(cols), rows), byrow = TRUE, nrow = rows)
mj <- matrix(rep(seq_len(cols), each = rows), byrow = TRUE, nrow = rows)
tmpl <- lapply(seq_len(ncirc), function(x) which((mi - centers[x, 1])^2 + (mj - centers[x, 2])^2 < radii[x]^2, arr.ind = TRUE))
inds <- do.call(rbind, tmpl)
grid <- matrix(1, nrow = rows, ncol = cols)
grid[inds] <- 2
alpha <- matrix(1e-5, ncol = cols, nrow = rows)
alpha[inds] <- 1e-7
## image(grid, asp=1)
## Define grid configuration for POET model
grid_setup <- list(
pqc_in_file = "./barite_fgcs_2.pqi",
pqc_db_file = "../barite/db_barite.dat", ## database file
grid_def = grid, ## grid definition, IDs according to the Phreeqc input
grid_size = c(dim_cols, dim_rows), ## grid size in meters
constant_cells = c() ## IDs of cells with constant concentration
)
bound_length <- cols / 10
bound_N <- list(
"type" = rep("constant", bound_length),
"sol_id" = rep(3, bound_length),
"cell" = seq(1, bound_length)
)
bound_W <- list(
"type" = rep("constant", bound_length),
"sol_id" = rep(3, bound_length),
"cell" = seq(1, bound_length)
)
bound_E <- list(
"type" = rep("constant", bound_length),
"sol_id" = rep(4, bound_length),
"cell" = seq(rows - bound_length + 1, rows)
)
bound_S <- list(
"type" = rep("constant", bound_length),
"sol_id" = rep(4, bound_length),
"cell" = seq(cols - bound_length + 1, cols)
)
diffusion_setup <- list(
boundaries = list(
"W" = bound_W,
"N" = bound_N,
"E" = bound_E,
"S" = bound_S
),
alpha_x = alpha,
alpha_y = alpha
)
dht_species <- c(
"H" = 7,
"O" = 7,
"Ba" = 7,
"Cl" = 7,
"S" = 7,
"Sr" = 7,
"Barite" = 4,
"Celestite" = 4
)
pht_species <- c(
"Ba" = 4,
"Cl" = 3,
"S" = 3,
"Sr" = 3,
"Barite" = 0,
"Celestite" = 0
)
chemistry_setup <- list(
dht_species = dht_species,
pht_species = pht_species
)
## 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
)

49
bench/fgcs/barite_fgcs_2.pqi
View File

@ -1,49 +0,0 @@
SOLUTION 1
units mol/kgw
water 1
temperature 25
pH 7.008
pe 10.798
S 6.205e-04
Sr 6.205e-04
END
SOLUTION 2
units mol/kgw
water 1
temperature 25
pH 7.008
pe 10.798
S 6.205e-04
Sr 6.205e-04
KINETICS 2
Barite
-m 0.00
-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
SOLUTION 4
units mol/kgw
water 1
temperature 25
Ba 0.2
Cl 0.4
END
RUN_CELLS
-cells 1 2 3 4
END

7
bench/fgcs/barite_fgcs_2_rt.R
View File

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

27
bench/surfex/CMakeLists.txt
View File

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

24
bench/surfex/ExBase.pqi
View File

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

40
bench/surfex/PoetEGU_surfex_500.R
View File

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

11
bench/surfex/PoetEGU_surfex_500_rt.R
View File

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

2
bench/surfex/README.org
View File

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

24
bench/surfex/SurfExBase.pqi
View File

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

108
bench/surfex/SurfexEGU.pqi
View File

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

157
bench/surfex/ex.R
View File

@ -1,37 +1,140 @@
rows <- 100
cols <- 100
## Time-stamp: "Last modified 2023-08-02 13:59:35 mluebke"
grid_def <- matrix(1, nrow = rows, ncol = cols)
database <- normalizePath("./SMILE_2021_11_01_TH.dat")
input_script <- normalizePath("./ExBase.pqi")
# 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
)
cat(paste(":: R This is a test 1\n"))
bound_def <- list(
"type" = rep("constant", cols),
"sol_id" = rep(2, cols),
"cell" = seq(1, cols)
)
#################################################################
## Section 1 ##
## Grid initialization ##
#################################################################
diffusion_setup <- list(
boundaries = list(
"N" = bound_def
),
alpha_x = 1e-6,
alpha_y = 1e-6
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"
)
chemistry_setup <- list()
##################################################################
## Section 2 ##
## Diffusion parameters and boundary conditions ##
##################################################################
vecinj_diffu <- list(
list(H = 0.147659686316291,
O = 0.0739242798146046,
Charge = 7.46361643222701e-20,
`C(-4)` = 2.92438561098248e-21,
`C(4)` = 2.65160558871092e-06,
Ca = 2.89001071336443e-05,
Cl = 0.000429291158114428,
`Fe(2)` = 1.90823391198114e-07,
`Fe(3)` = 3.10832423034763e-12,
`H(0)` = 2.7888235127385e-15,
K = 2.5301787e-06,
Mg = 2.31391999937907e-05,
Na = 0.00036746969,
`S(-2)` = 1.01376078438546e-14,
`S(2)` = 1.42247026981542e-19,
`S(4)` = 9.49422092568557e-18,
`S(6)` = 2.19812504654191e-05,
Sr = 6.01218519999999e-07,
`U(4)` = 4.82255946569383e-12,
`U(5)` = 5.49050615347901e-13,
`U(6)` = 1.32462838991902e-09)
)
vecinj <- do.call(rbind.data.frame, vecinj_diffu)
names(vecinj) <- grid$props
## diffusion coefficients
alpha_diffu <- c(H = 1E-6, O = 1E-6, Charge = 1E-6, `C(-4)` = 1E-6,
`C(4)` = 1E-6, Ca = 1E-6, Cl = 1E-6, `Fe(2)` = 1E-6,
`Fe(3)` = 1E-6, `H(0)` = 1E-6, K = 1E-6, Mg = 1E-6,
Na = 1E-6, `S(-2)` = 1E-6, `S(2)` = 1E-6,
`S(4)` = 1E-6, `S(6)` = 1E-6, Sr = 1E-6,
`U(4)` = 1E-6, `U(5)` = 1E-6, `U(6)` = 1E-6)
## list of boundary conditions/inner nodes
## vecinj_inner <- list(
## list(1,1,1)
## )
boundary <- list(
"N" = rep(1, n),
"E" = rep(0, n),
"S" = rep(0, n),
"W" = rep(0, n)
)
diffu_list <- names(alpha_diffu)
vecinj <- do.call(rbind.data.frame, vecinj_diffu)
names(vecinj) <- names(init_cell)
diffusion <- list(
init = as.data.frame(init_cell, check.names = FALSE),
vecinj = vecinj,
# vecinj_inner = vecinj_inner,
vecinj_index = boundary,
alpha = alpha_diffu
)
#################################################################
## Section 3 ##
## Chemistry module (Phreeqc) ##
#################################################################
chemistry <- list(
database = database,
input_script = input_script
)
#################################################################
## Section 4 ##
## Putting all those things together ##
#################################################################
iterations <- 10
dt <- 200
# Define a setup list for simulation configuration
setup <- list(
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
grid = grid,
diffusion = diffusion,
chemistry = chemistry,
iterations = iterations,
timesteps = rep(dt, iterations),
store_result = TRUE
)

7
bench/surfex/ex_rt.R
View File

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

158
bench/surfex/surfex.R
View File

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

10
bench/surfex/surfex_rt.R
View File

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

8
bin/RunDaosBenchmarks.sh Normal file
View File

@ -0,0 +1,8 @@
taskset --cpu-list 0-23:2 mpirun --allow-run-as-root -x DAOS_POOL=test_pool -n 12 ./poet ../share/poet/bench/daos/dolo_diffu_edge.R ../../POETR >> ../../POET.out
# taskset --cpu-list 0-23:2 mpirun --allow-run-as-root -x DAOS_POOL=test_pool -n $i ./build/src/kivibench-DAOSKV -x 10000 -y 3 -m 1000 -k 10 -v 12000 --csv >> benchmarks/clientbig$i.csv

BIN
bin/poet Executable file
View File

Binary file not shown.

572
docs/20221216_Scheme_PORT_en.svg Normal file
View File

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

After

Width:  |  Height:  |  Size: 73 KiB

568
docs/20230720_Scheme_POET_en.svg Normal file
View File

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

After

Width:  |  Height:  |  Size: 73 KiB

3
docs/CMakeLists.txt
View File

@ -12,8 +12,9 @@ if(DOXYGEN_FOUND)
set(DOXYGEN_PROJECT_NUMBER ${POET_VERSION})
doxygen_add_docs(doxygen
${PROJECT_SOURCE_DIR}/src
${PROJECT_SOURCE_DIR}/include
${PROJECT_SOURCE_DIR}/README.md
${PROJECT_SOURCE_DIR}/docs/Input_Scripts.md
${PROJECT_SOURCE_DIR}/docs/Output.md
COMMENT "Generate html pages")
endif()

86
docs/Input_Scripts.md Normal file
View File

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

50
docs/Output.md
View File

@ -36,49 +36,33 @@ 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 |
| chemistry | measured time in *chemistry* subroutine |
| diffusion | measured time in *diffusion* subroutine |
| simtime\_transport | measured time in *transport* subroutine |
| simtime\_chemistry | measured time in *chemistry* subroutine (actual parallelized part) |
### 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 |
|-----------------------|-----------------------------------------------------------|
| 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\_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 |
| dht\_get\_time | time to retreive data from DHT |
| dh\_hits | count of data points retrieved from DHT |
| dht\_miss | count of misses/count of data points written to DHT |
| dht\_evictions | count of data points evicted by another write operation |

281
docs/POET.drawio
View File

File diff suppressed because one or more lines are too long

4
docs/POET_scheme.svg
View File

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

Before

Width:  |  Height:  |  Size: 813 KiB

602
docs/Scheme_POET_en.svg
View File

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

Before

Width:  |  Height:  |  Size: 127 KiB

1
ext/doctest Submodule

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

1
ext/litephreeqc

@ -1 +0,0 @@
Subproject commit 953c752431d2b2758268083f407f943843efc7ad

1
ext/phreeqcrm Submodule

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

2
ext/tug

@ -1 +1 @@
Subproject commit 9c4aeee410c71d064f7567143d4f8d6451ade75a
Subproject commit 25855da6b2930559b542bbadb16299932332d6a3

152
src/Chemistry/ChemistryModule.hpp → include/poet/ChemistryModule.hpp
View File

@ -1,19 +1,17 @@
// Time-stamp: "Last modified 2023-08-15 14:36:28 mluebke"
#ifndef CHEMISTRYMODULE_H_
#define CHEMISTRYMODULE_H_
#include "DataStructures/Field.hpp"
#include "DataStructures/NamedVector.hpp"
#include "DHT_Wrapper.hpp"
#include "DataStructures.hpp"
#include "Interpolation.hpp"
#include "IrmResult.h"
#include "PhreeqcRM.h"
#include "SimParams.hpp"
#include "ChemistryDefs.hpp"
#include "Init/InitialList.hpp"
#include "NameDouble.h"
#include "SurrogateModels/DHT_Wrapper.hpp"
#include "SurrogateModels/Interpolation.hpp"
#include "PhreeqcRunner.hpp"
#include <array>
#include <cstddef>
#include <cstdint>
#include <map>
#include <memory>
@ -26,7 +24,7 @@ namespace poet {
* \brief Wrapper around PhreeqcRM to provide POET specific parallelization with
* easy access.
*/
class ChemistryModule {
class ChemistryModule : public PhreeqcRM {
public:
/**
* Creates a new instance of Chemistry module with given grid cell count, work
@ -43,80 +41,61 @@ 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 wp_size,
const InitialList::ChemistryInit chem_params,
MPI_Comm communicator);
ChemistryModule(uint32_t nxyz, uint32_t wp_size, std::uint32_t maxiter,
const ChemistryParams &chem_param, MPI_Comm communicator);
/**
* Deconstructor, which frees DHT data structure if used.
*/
~ChemistryModule();
void masterSetField(Field field);
/**
* Parses the input script and extract information needed during runtime.
*
* **Only run by master**.
*
* Database must be loaded beforehand.
*
* \param input_script_path Path to input script to parse.
*/
void RunInitFile(const std::string &input_script_path);
/**
* Run the chemical simulation with parameters set.
*/
void simulate(double dt);
void RunCells();
/**
* 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; }
void setFilePadding(std::uint32_t maxiter) {
this->file_pad =
static_cast<std::uint8_t>(std::ceil(std::log10(maxiter + 1)));
}
/**
* 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);
struct SurrogateSetup {
std::vector<std::string> prop_names;
std::array<double, 2> base_totals;
bool has_het_ids;
bool dht_enabled;
std::uint32_t dht_size_mb;
int dht_snaps;
std::string dht_out_dir;
bool interp_enabled;
std::uint32_t interp_bucket_size;
std::uint32_t interp_size_mb;
std::uint32_t interp_min_entries;
bool ai_surrogate_enabled;
};
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;
this->ai_surrogate_enabled = setup.ai_surrogate_enabled;
this->base_totals = setup.base_totals;
if (this->dht_enabled || this->interp_enabled) {
this->initializeDHT(setup.dht_size_mb, this->params.dht_species,
setup.has_het_ids);
if (setup.dht_snaps != DHT_SNAPS_DISABLED) {
this->setDHTSnapshots(setup.dht_snaps, setup.dht_out_dir);
}
}
if (this->interp_enabled) {
this->initializeInterp(setup.interp_bucket_size, setup.interp_size_mb,
setup.interp_min_entries,
this->params.interp_species);
}
}
/**
* Default work package size.
*/
static constexpr uint32_t CHEM_DEFAULT_WP_SIZE = 5;
/**
* Intended to alias input parameters for grid initialization with a single
@ -139,6 +118,17 @@ 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.
*/
@ -233,11 +223,6 @@ public:
this->print_progessbar = enabled;
};
/**
* **Master only** Set the ai surrogate validity vector from R
*/
void set_ai_surrogate_validity_vector(std::vector<int> r_vector);
std::vector<uint32_t> GetWorkerInterpolationCalls() const;
std::vector<double> GetWorkerInterpolationWriteTimings() const;
@ -247,12 +232,9 @@ public:
std::vector<uint32_t> GetWorkerPHTCacheHits() const;
std::vector<int> ai_surrogate_validity_vector;
protected:
void initializeDHT(uint32_t size_mb,
const NamedVector<std::uint32_t> &key_species,
bool has_het_ids);
const NamedVector<std::uint32_t> &key_species);
void setDHTSnapshots(int type, const std::string &out_dir);
void setDHTReadFile(const std::string &input_file);
@ -261,9 +243,12 @@ protected:
const NamedVector<std::uint32_t> &key_species);
enum {
CHEM_FIELD_INIT,
CHEM_INIT,
CHEM_WP_SIZE,
CHEM_INIT_SPECIES,
CHEM_DHT_ENABLE,
CHEM_DHT_SIGNIF_VEC,
CHEM_DHT_PROP_TYPE_VEC,
CHEM_DHT_SNAPS,
CHEM_DHT_READ_FILE,
CHEM_IP_ENABLE,
@ -271,8 +256,7 @@ protected:
CHEM_IP_SIGNIF_VEC,
CHEM_WORK_LOOP,
CHEM_PERF,
CHEM_BREAK_MAIN_LOOP,
CHEM_AI_BCAST_VALIDITY
CHEM_BREAK_MAIN_LOOP
};
enum { LOOP_WORK, LOOP_END };
@ -311,7 +295,7 @@ protected:
using worker_list_t = std::vector<struct worker_info_s>;
using workpointer_t = std::vector<double>::iterator;
void MasterRunParallel(double dt);
void MasterRunParallel();
void MasterRunSequential();
void MasterSendPkgs(worker_list_t &w_list, workpointer_t &work_pointer,
@ -337,7 +321,7 @@ protected:
void WorkerPerfToMaster(int type, const struct worker_s &timings);
void WorkerMetricsToMaster(int type);
void WorkerRunWorkPackage(WorkPackage &work_package, double dSimTime,
IRM_RESULT WorkerRunWorkPackage(WorkPackage &work_package, double dSimTime,
double dTimestep);
std::vector<uint32_t> CalculateWPSizesVector(uint32_t n_cells,
@ -361,7 +345,7 @@ protected:
bool is_sequential;
bool is_master;
uint32_t wp_size;
uint32_t wp_size{CHEM_DEFAULT_WP_SIZE};
bool dht_enabled{false};
int dht_snaps_type{DHT_SNAPS_DISABLED};
std::string dht_file_out_dir;
@ -371,8 +355,6 @@ protected:
bool interp_enabled{false};
std::unique_ptr<poet::InterpolationModule> interp;
bool ai_surrogate_enabled{false};
static constexpr uint32_t BUFFER_OFFSET = 5;
inline void ChemBCast(void *buf, int count, MPI_Datatype datatype) const {
@ -391,7 +373,7 @@ protected:
bool print_progessbar{false};
std::uint8_t file_pad{1};
std::uint32_t file_pad;
double chem_t{0.};
@ -401,9 +383,11 @@ protected:
Field chem_field;
const InitialList::ChemistryInit params;
static constexpr int MODULE_COUNT = 5;
std::unique_ptr<PhreeqcRunner> pqc_runner;
const ChemistryParams &params;
std::array<std::uint32_t, MODULE_COUNT> speciesPerModule{};
};
} // namespace poet

46
src/Chemistry/SurrogateModels/DHT.h → include/poet/DHT.h
View File

@ -117,9 +117,6 @@ typedef struct {
unsigned int index_count;
int (*accumulate_callback)(int, void *, int, void *);
size_t sum_idx;
size_t cnt_idx;
#ifdef DHT_STATISTICS
/** Detailed statistics of the usage of the DHT. */
DHT_stats *stats;
@ -131,8 +128,7 @@ extern void DHT_set_accumulate_callback(DHT *table,
void *));
extern int DHT_write_accumulate(DHT *table, const void *key, int send_size,
void *data, uint32_t *proc, uint32_t *index,
int *callback_ret);
void *data, uint32_t *proc, uint32_t *index, int *callback_ret);
/**
* @brief Create a DHT.
@ -288,8 +284,44 @@ extern int DHT_free(DHT *table, int *eviction_counter, int *readerror_counter);
*/
extern int DHT_print_statistics(DHT *table);
extern float DHT_get_used_idx_factor(DHT *table, int with_reset);
/**
* @brief Determine destination rank and index.
*
* This is done by looping over all possbile indices. First of all, set a
* temporary index to zero and copy count of bytes for each index into the
* memory area of the temporary index. After that the current index is
* calculated by the temporary index modulo the table size. The destination rank
* of the process is simply determined by hash modulo the communicator size.
*
* @param hash Calculated 64 bit hash.
* @param comm_size Communicator size.
* @param table_size Count of buckets per process.
* @param dest_rank Reference to the destination rank variable.
* @param index Pointer to the array index.
* @param index_count Count of possible indeces.
*/
static void determine_dest(uint64_t hash, int comm_size,
unsigned int table_size, unsigned int *dest_rank,
uint64_t *index, unsigned int index_count);
extern int DHT_flush(DHT *table);
/**
* @brief Set the occupied flag.
*
* This will set the first bit of a bucket to 1.
*
* @param flag_byte First byte of a bucket.
*/
static void set_flag(char *flag_byte);
/**
* @brief Get the occupied flag.
*
* This function determines whether the occupied flag of a bucket was set or
* not.
*
* @param flag_byte First byte of a bucket.
* @return int Returns 1 for true or 0 for false.
*/
static int read_flag(char flag_byte);
#endif /* DHT_H */

23
src/Chemistry/SurrogateModels/DHT_Wrapper.hpp → include/poet/DHT_Wrapper.hpp
View File

@ -1,3 +1,4 @@
// Time-stamp: "Last modified 2023-09-08 14:43:02 mluebke"
/*
** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
@ -22,18 +23,19 @@
#ifndef DHT_WRAPPER_H
#define DHT_WRAPPER_H
#include "Base/RInsidePOET.hpp"
#include "DataStructures/NamedVector.hpp"
#include "Chemistry/ChemistryDefs.hpp"
#include "Init/InitialList.hpp"
#include "DataStructures.hpp"
#include "LookupKey.hpp"
#include "RInsidePOET.hpp"
#include "SimParams.hpp"
#include "enums.hpp"
#include "poet/HashFunctions.hpp"
#include "poet/LookupKey.hpp"
#include "poet/Rounding.hpp"
#include <array>
#include <cstdint>
#include <limits>
#include <string>
#include <unordered_map>
#include <utility>
#include <vector>
@ -86,8 +88,8 @@ public:
const NamedVector<std::uint32_t> &key_species,
const std::vector<std::int32_t> &key_indices,
const std::vector<std::string> &output_names,
const InitialList::ChemistryHookFunctions &hooks,
uint32_t data_count, bool with_interp, bool has_het_ids);
const ChemistryParams::Chem_Hook_Functions &hooks,
uint32_t data_count, bool with_interp);
/**
* @brief Destroy the dht wrapper object
*
@ -258,13 +260,12 @@ private:
const std::vector<std::string> &output_names;
const InitialList::ChemistryHookFunctions &hooks;
const ChemistryParams::Chem_Hook_Functions &hooks;
const bool with_interp;
DHT_ResultObject dht_results;
std::array<double, 2> base_totals{0};
bool has_het_ids{false};
};
} // namespace poet

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

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

42
src/DataStructures/Field.hpp → include/poet/DataStructures.hpp
View File

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

52
src/Transport/DiffusionModule.hpp → include/poet/DiffusionModule.hpp
View File

@ -21,13 +21,19 @@
#ifndef DIFFUSION_MODULE_H
#define DIFFUSION_MODULE_H
#include "DataStructures/Field.hpp"
#include "Init/InitialList.hpp"
#include <sys/types.h>
#include "DataStructures.hpp"
#include "SimParams.hpp"
#include "poet/SimParams.hpp"
#include <array>
#include <cmath>
#include <cstdint>
#include <poet/Grid.hpp>
#include <string>
#include <tug/BoundaryCondition.hpp>
#include <tug/Diffusion.hpp>
#include <vector>
namespace poet {
/**
* @brief Class describing transport simulation
*
@ -46,8 +52,8 @@ public:
*
* @param R RRuntime object
*/
DiffusionModule(const InitialList::DiffusionInit &init_list, Field field)
: param_list(init_list), transport_field(field){};
DiffusionModule(const poet::DiffusionParams &diffu_args,
const poet::GridParams &grid_params);
/**
* @brief Run simulation for one iteration
@ -57,6 +63,14 @@ public:
*/
void simulate(double dt);
/**
* @brief End simulation
*
* All measured timings are distributed to the R runtime
*
*/
void end();
/**
* @brief Get the transport time
*
@ -69,23 +83,39 @@ public:
*
* \return Reference to the diffusion field.
*/
Field &getField() { return this->transport_field; }
Field &getField() { return this->t_field; }
private:
/**
* @brief Instance of RRuntime
*
*/
// RRuntime &R;
InitialList::DiffusionInit param_list;
enum { DIM_1D = 1, DIM_2D };
Field transport_field;
void initialize(const poet::DiffusionParams &args,
std::uint32_t n_grid_cells);
uint8_t dim;
uint32_t prop_count;
tug::diffusion::TugInput diff_input;
std::vector<double> alpha;
std::vector<uint32_t> index_constant_cells;
std::vector<std::string> prop_names;
std::vector<tug::bc::BoundaryCondition> bc_vec;
Field t_field;
uint32_t n_cells_per_prop;
/**
* @brief time spent for transport
*
*/
double transport_t = 0.;
double transport_t = 0.f;
};
} // namespace poet

116
include/poet/Grid.hpp Normal file
View File

@ -0,0 +1,116 @@
/*
** 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 "poet/SimParams.hpp"
#include <RInside.h>
#include <Rcpp.h>
#include <array>
#include <cstddef>
#include <cstdint>
#include <list>
#include <map>
#include <string>
#include <vector>
#define MAX_DIM 2
namespace poet {
enum { GRID_X_DIR, GRID_Y_DIR, GRID_Z_DIR };
using StateMemory = struct s_StateMemory {
std::vector<double> mem;
std::vector<std::string> props;
};
using FlowInputOutputInfo = struct inOut_info {
std::string input_field;
std::string output_field;
};
constexpr const char *GRID_MODULE_NAME = "grid_init";
/**
* @brief Class describing the grid
*
* Providing methods to shuffle and unshuffle grid (for the master) as also to
* import and export a work package (for worker).
*
* @todo find better abstraction
*
*/
class Grid {
public:
Grid();
~Grid();
void InitModuleFromParams(const poet::GridParams &grid_args);
void SetGridDimension(uint8_t dim);
void SetGridCellCount(uint32_t n_x, uint32_t n_y = 0, uint32_t n_z = 0);
void SetGridSize(double s_x, double s_y = 0., double s_z = 0.);
void SetPropNames(const std::vector<std::string> &prop_names);
void PushbackModuleFlow(const std::string &input, const std::string &output);
void PreModuleFieldCopy(uint32_t tick);
void InitGridFromScratch(const Rcpp::DataFrame &init_cell);
void InitGridFromRDS();
void InitGridFromPhreeqc();
auto GetGridDimension() const -> uint8_t;
auto GetTotalCellCount() const -> uint32_t;
auto GetGridCellsCount(uint8_t direction) const -> uint32_t;
auto GetGridSize(uint8_t direction) const -> uint32_t;
auto RegisterState(std::string module_name, std::vector<std::string> props)
-> StateMemory *;
auto GetStatePointer(std::string module_name) -> StateMemory *;
auto GetInitialGrid() const -> StateMemory *;
auto GetSpeciesCount() const -> uint32_t;
auto GetPropNames() const -> std::vector<std::string>;
auto GetSpeciesByName(std::string name,
std::string module_name = poet::GRID_MODULE_NAME) const
-> std::vector<double>;
void WriteFieldsToR(RInside &R);
private:
std::vector<FlowInputOutputInfo> flow_vec;
std::uint8_t dim = 0;
std::array<double, MAX_DIM> grid_size;
std::array<std::uint32_t, MAX_DIM> n_cells;
std::map<std::string, StateMemory *> state_register;
StateMemory *grid_init = std::nullptr_t();
std::vector<std::string> prop_names;
};
} // namespace poet
#endif // GRID_H

1
src/Chemistry/SurrogateModels/HashFunctions.hpp → include/poet/HashFunctions.hpp
View File

@ -1,3 +1,4 @@
// // Time-stamp: "Last modified 2023-03-31 14:59:49 mluebke"
/*
** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of

26
src/Chemistry/SurrogateModels/Interpolation.hpp → include/poet/Interpolation.hpp
View File

@ -1,25 +1,29 @@
// Time-stamp: "Last modified 2023-08-16 16:49:31 mluebke"
#ifndef INTERPOLATION_H_
#define INTERPOLATION_H_
#include "DataStructures/NamedVector.hpp"
#include "DHT.h"
#include "DHT_Wrapper.hpp"
#include "Init/InitialList.hpp"
#include "DataStructures.hpp"
#include "LookupKey.hpp"
#include "Rounding.hpp"
#include "poet/DHT_Wrapper.hpp"
#include "poet/Rounding.hpp"
#include "poet/SimParams.hpp"
#include <cassert>
#include <iostream>
#include <list>
#include <memory>
#include <mpi.h>
#include <string>
#include <utility>
extern "C" {
#include "DHT.h"
#include "poet/DHT.h"
}
#include "poet/LookupKey.hpp"
#include <cstdint>
#include <functional>
#include <unordered_map>
#include <vector>
@ -162,12 +166,10 @@ public:
const NamedVector<std::uint32_t> &interp_key_signifs,
const std::vector<std::int32_t> &dht_key_indices,
const std::vector<std::string> &out_names,
const InitialList::ChemistryHookFunctions &hooks);
const ChemistryParams::Chem_Hook_Functions &hooks);
enum result_status { RES_OK, INSUFFICIENT_DATA, NOT_NEEDED };
DHT *getDHTObject() { return this->pht->getDHTObject(); }
struct InterpolationResult {
std::vector<std::vector<double>> results;
std::vector<result_status> status;
@ -260,11 +262,9 @@ private:
return out_key;
}
const InitialList::ChemistryHookFunctions &hooks;
const ChemistryParams::Chem_Hook_Functions &hooks;
const std::vector<std::string> &out_names;
const std::vector<std::string> dht_names;
std::unordered_map<int, std::vector<std::int32_t>> to_calc_cache;
};
} // namespace poet

20
src/Chemistry/SurrogateModels/LookupKey.hpp → include/poet/LookupKey.hpp
View File

@ -1,30 +1,18 @@
// Time-stamp: "Last modified 2023-08-11 10:12:52 mluebke"
#ifndef LOOKUPKEY_H_
#define LOOKUPKEY_H_
#include "HashFunctions.hpp"
#include "poet/HashFunctions.hpp"
#include <cstdint>
#include <cstring>
#include <vector>
namespace poet {
constexpr std::int8_t SC_NOTATION_EXPONENT_MASK = -128;
constexpr std::int64_t SC_NOTATION_SIGNIFICANT_MASK = 0xFFFFFFFFFFFF;
struct Lookup_SC_notation {
std::int8_t exp : 8;
std::int64_t significant : 56;
constexpr static Lookup_SC_notation nan() {
return {SC_NOTATION_EXPONENT_MASK, SC_NOTATION_SIGNIFICANT_MASK};
}
constexpr bool isnan() const {
return !!(exp == SC_NOTATION_EXPONENT_MASK &&
significant == SC_NOTATION_SIGNIFICANT_MASK);
}
};
union Lookup_Keyelement {
@ -35,10 +23,6 @@ union Lookup_Keyelement {
return std::memcmp(this, &other, sizeof(Lookup_Keyelement)) == 0 ? true
: false;
}
template <typename T> bool operator>(const T &other) const {
return this->sc_notation.significant > other;
}
};
class LookupKey : public std::vector<Lookup_Keyelement> {

1
src/Base/Macros.hpp → include/poet/Macros.hpp
View File

@ -1,3 +1,4 @@
// Time-stamp: "Last modified 2023-08-09 14:16:04 mluebke"
#ifndef MACROS_H
#define MACROS_H

59
include/poet/RInsidePOET.hpp Normal file
View File

@ -0,0 +1,59 @@
#ifndef RPOET_H_
#define RPOET_H_
#include <RInside.h>
#include <Rcpp.h>
#include <cstddef>
#include <exception>
#include <optional>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
class RInsidePOET : public RInside {
public:
static RInsidePOET &getInstance() {
static RInsidePOET instance;
return instance;
}
RInsidePOET(RInsidePOET const &) = delete;
void operator=(RInsidePOET const &) = delete;
inline bool checkIfExists(const std::string &R_name,
const std::string &where) {
return Rcpp::as<bool>(
this->parseEval("'" + R_name + "' %in% names(" + where + ")"));
}
private:
RInsidePOET() : RInside(){};
};
template <typename T> class RHookFunction {
public:
RHookFunction() {}
RHookFunction(RInside &R, const std::string &f_name) {
try {
this->func = Rcpp::Function(Rcpp::as<SEXP>(R.parseEval(f_name.c_str())));
} catch (const std::exception &e) {
}
}
template <typename... Args> T operator()(Args... args) const {
if (func.has_value()) {
return (Rcpp::as<T>(this->func.value()(args...)));
} else {
throw std::exception();
}
}
bool isValid() const { return this->func.has_value(); }
private:
std::optional<Rcpp::Function> func;
};
#endif // RPOET_H_

13
src/Chemistry/SurrogateModels/Rounding.hpp → include/poet/Rounding.hpp
View File

@ -20,11 +20,6 @@ class DHT_Rounder {
public:
Lookup_Keyelement round(const double &value, std::uint32_t signif,
bool is_ho) {
if (std::isnan(value)) {
return {.sc_notation = Lookup_SC_notation::nan()};
}
std::int8_t exp =
static_cast<std::int8_t>(std::floor(std::log10(std::fabs(value))));
@ -65,14 +60,6 @@ public:
std::uint32_t signif) {
Lookup_Keyelement new_val = value;
if (value.sc_notation.isnan()) {
return {.sc_notation = Lookup_SC_notation::nan()};
}
if (signif == 0) {
return {.sc_notation = {0, value > 0}};
}
std::uint32_t diff_signif =
static_cast<std::uint32_t>(
std::ceil(std::log10(std::abs(value.sc_notation.significant)))) -

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