max/poet
1
0
Fork 0
mirror of https://git.gfz-potsdam.de/naaice/poet.git synced 2025-12-16 04:48:23 +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:v0.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

No commits in common. "main" and "v0.2" have entirely different histories.
main ... v0.2

139 changed files with 5748 additions and 9414 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"
}

4
.gitignore vendored
View File

@ -140,7 +140,3 @@ vignettes/*.pdf
# End of https://www.toptal.com/developers/gitignore/api/c,c++,r,cmake
build/
/.cache/
.vscode
.codechecker

151
.gitlab-ci.yml
View File

@ -16,104 +16,99 @@
# This specific template is located at:
# https://gitlab.com/gitlab-org/gitlab/-/blob/master/lib/gitlab/ci/templates/Getting-Started.gitlab-ci.yml
image: git.gfz-potsdam.de:5000/naaice/poet:ci
image: git.gfz-potsdam.de:5000/sec34/port:builder
stages: # List of stages for jobs, and their order of execution
- build
- release
- test
variables:
GIT_SUBMODULE_STRATEGY: recursive
SOURCE_ARCHIVE_NAME: 'poet_${CI_COMMIT_TAG}_sources.tar.gz'
SOURCE_ARCHIVE_NAME: 'port_${CI_COMMIT_TAG}_sources.tar.gz'
CHANGELOG_FILE: 'commit_changelog.md'
test: # This job runs in the build stage, which runs first.
build-poet: # This job runs in the build stage, which runs first.
stage: build
script:
- mkdir build && cd build
- cmake ..
- make -j$(nproc)
rules:
- if: $CI_PIPELINE_SOURCE == 'merge_request_event'
test-poet:
stage: test
script:
- mkdir -p build && cd build
- cmake -DPOET_ENABLE_TESTING=ON -DPOET_PREPROCESS_BENCHS=OFF -DCMAKE_BUILD_TYPE=Release ..
- mkdir build_test && cd build_test
- cmake -DPOET_ENABLE_TESTING=ON ..
- make -j$(nproc) check
rules:
- if: $CI_PIPELINE_SOURCE == 'merge_request_event'
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}
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: 'PORT $CI_COMMIT_TAG'
description: ${CHANGELOG_FILE}
assets:
links:
- name: '${SOURCE_ARCHIVE_NAME}'
url: 'https://git.gfz-potsdam.de/sec34/port/-/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 {public,build}
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

11
.gitmodules vendored
View File

@ -1,7 +1,10 @@
[submodule "ext/tug"]
path = ext/tug
url = ../tug.git
url = ../../sec34/tug.git
[submodule "ext/litephreeqc"]
path = ext/litephreeqc
url = ../litephreeqc.git
[submodule "ext/phreeqcrm"]
path = ext/phreeqcrm
url = ../../mluebke/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

52
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(
COMMAND ${R_EXE} RHOME
OUTPUT_VARIABLE R_ROOT_DIR
OUTPUT_STRIP_TRAILING_WHITESPACE
)
if(R_EXE)
execute_process(
COMMAND ${R_EXE} RHOME
OUTPUT_VARIABLE R_ROOT_DIR
OUTPUT_STRIP_TRAILING_WHITESPACE
)
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
)
find_path(
R_INCLUDE_DIR R.h
HINTS ${R_ROOT_DIR}
PATHS /usr/include /usr/local/include /usr/share
PATH_SUFFIXES include/R R/include
)
find_library(
R_LIBRARY libR.so
HINTS ${R_ROOT_DIR}/lib
REQUIRED
)
find_library(
R_LIBRARY libR.so
HINTS ${R_ROOT_DIR}/lib
)
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)

6
CMake/POET_Scripts.cmake
View File

@ -9,11 +9,11 @@ macro(get_POET_version)
OUTPUT_VARIABLE POET_GIT_BRANCH
OUTPUT_STRIP_TRAILING_WHITESPACE)
execute_process(
COMMAND ${GIT_EXECUTABLE} describe --tags --always
COMMAND ${GIT_EXECUTABLE} describe --always
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}")

30
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,28 @@ 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(data)
add_subdirectory(app)
add_subdirectory(bench/dolo_diffu_inner)
# 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.

334
README.md
View File

@ -1,108 +1,81 @@
<!--
Time-stamp: "Last modified 2023-01-19 12:06:10 delucia"
-->
# Forked Project
*PORT* is a fork of [POET](https://doi.org/10.5281/zenodo.4757913)
integrating a standalone component for transport computations and
leveraging PHREEQC_RM as geochemical solver. The following README is
also applicable for this project.
![PORT's Coupling Scheme](./docs/20221216_Scheme_PORT_en.svg)
# 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's Coupling Scheme](./docs/POET_scheme.svg)
POET is a coupled reactive transport simulator implementing a parallel
architecture and a fast, original MPI-based Distributed Hash Table.
## 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://sec34.git-pages.gfz-potsdam.de/port).
## 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_USE_PRM_BACKEND**=_bollean_ - use the PhreeqcRM parallelization instead
of POET's one. Intended for debugging purposes for modellers.
### Example: Build from scratch
@ -115,7 +88,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 +96,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,64 +111,56 @@ The correspondending directory tree would look like this:
```sh
poet
├── bin
│   ├── poet
│   └── poet_init
│ └── poet
├── R_lib
│ └── kin_r_library.R
└── share
└── poet
├── barite
│   ├── barite_200.qs2
│   ├── barite_200_rt.R
│   ├── barite_het.qs2
│   └── barite_het_rt.R
├── dolo
│   ├── dolo_inner_large.qs2
│   ├── dolo_inner_large_rt.R
│   ├── dolo_interp.qs2
│   └── dolo_interp_rt.R
└── surfex
├── PoetEGU_surfex_500.qs2
└── PoetEGU_surfex_500_rt.R
├── bench
│ ├── dolo_diffu_inner_large.R
│ ├── dolo_diffu_inner.R
│ └── dolo_inner.pqi
└── examples
├── dol.pqi
├── phreeqc_kin.dat
├── SimDol1D_diffu.R
└── SimDol2D_diffu.R
```
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
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 R libraries will be loaded at runtime and the paths are hardcoded
absolute paths inside `poet.cpp`. So, if you consider to move
`bin/poet` either change paths of the R source files and recompile
POET or also move `R_lib/*` relative to the binary.
## Running
Run POET by `mpirun ./poet [OPTIONS] <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/examples/SimDol2D_diffu.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_) |
| **--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 |
| **--dht-file=** | `<SNAPSHOT>` | initializes DHT with the given snapshot file |
| **--interp-size** | _1-n_ | size of PHT (interpolation) per process in megabyte |
| **--interp-bucket-entries** | _1-n_ | number of entries to store at maximum in one PHT bucket |
| **--interp-min** | _1-n_ | number of entries in PHT bucket needed to start interpolation |
| Option | Value | Description |
|--------------------------|--------------|--------------------------------------------------------------------------------------------------------------------------|
| **--work-package-size=** | _1..n_ | size of work packages (defaults to _5_) |
| **--ignore-result** | | disables store of simulation resuls |
| **--dht** | | enabling DHT usage (defaults to _OFF_) |
| **--dht-signif=** | _1..n_ | set rounding to number of significant digits (defaults to _5_) (it is recommended to use `signif_vec` in R input script) |
| **--dht-strategy=** | _0-1_ | change DHT strategy. **NOT IMPLEMENTED YET** (Defaults to _0_) |
| **--dht-size=** | _1-n_ | size of DHT per process involved in megabyte (defaults to _1000 MByte_) |
| **--dht-snaps=** | _0-2_ | disable or enable storage of DHT snapshots |
| **--dht-file=** | `<SNAPSHOT>` | initializes DHT with the given snapshot file |
#### Additions to `dht-signif`
Only used if no vector is given in setup file. For individual values
per column use R vector `signif_vector` in `SIMFILE`.
#### Additions to `dht-snaps`
@ -210,106 +175,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:
`SimDol2D_diffu.R`. As transport a simple fixed-coefficient diffusion is used.
It's a 2D, 100x100 grid, simulating 10 time steps. To start the simulation with
4 processes `cd` into your previously installed POET-dir
`<POET_INSTALL_DIR>/bin` and run:
```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/examples/SimDol2D_diffu.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/examples/SimDol2D_diffu.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 +205,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)
}

411
R_lib/kin_r_library.R
View File

@ -1,4 +1,4 @@
### Copyright (C) 2018-2023 Marco De Lucia, Max Luebke (GFZ 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
@ -13,199 +13,270 @@
### 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) {
## Setup the directory where we will store the results
if (!dir.exists(out_dir)) {
dir.create(out_dir)
msgm("created directory ", out_dir)
## 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
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)) {
msgm("store_result doesn't exist!")
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$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")
setup$iter <- 1
setup$maxiter <- setup$iterations
setup$timesteps <- setup$timesteps
setup$simulation_time <- 0
if (is.null(setup[["store_result"]])) {
setup$store_result <- TRUE
if (is.null(setup[["store_result"]])) {
setup$store_result <- TRUE
}
if (setup$store_result) {
if (is.null(setup[["out_save"]])) {
setup$out_save <- seq(1, setup$iterations)
}
}
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)
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) {
iter <- setup$iter
# print(iter)
## max digits for iterations
dgts <- as.integer(ceiling(log10(setup$maxiter + 1)))
## string format to use in sprintf
fmt <- paste0("%0", dgts, "d")
## Write on disk state_T and state_C after every iteration
## comprised in setup$out_save
if (setup$store_result) {
if (iter %in% setup$out_save) {
nameout <- paste0(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
)
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)
msgm("results stored in <", nameout, ">")
}
master_iteration_end <- function(setup) {
iter <- setup$iter
## MDL Write on disk state_T and state_C after every iteration
## comprised in setup$out_save
if (setup$store_result) {
if (iter %in% setup$out_save) {
nameout <- paste0(fileout, "/iter_", sprintf("%03d", iter), ".rds")
info <- list(
tr_req_dt = as.integer(setup$req_dt)
# tr_allow_dt = setup$allowed_dt,
# tr_inniter = as.integer(setup$inniter)
)
saveRDS(list(
T = setup$state_T, C = setup$state_C,
simtime = as.integer(setup$simtime),
tr_info = info
), file = nameout)
msgm("results stored in <", nameout, ">")
}
## Add last time step to simulation time
setup$simulation_time <- setup$simulation_time + setup$timesteps[iter]
## msgm("done iteration", iter, "/", length(setup$timesteps))
setup$iter <- setup$iter + 1
return(setup)
}
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()
}
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 = 3)
# names(to_store) <- c("Sim", "Flow", "Transport", "Chemistry", "DHT")
names(to_store) <- c("Sim", "Transport", "DHT")
## 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
}
saveRDS(to_store, file = paste0(fileout, "/setup.rds"))
msgm("initialization stored in ", paste0(fileout, "/setup.rds"))
}
GetWorkPackageSizesVector <- function(n_packages, package_size, len) {
ids <- rep(1:n_packages, times = package_size, each = 1)[1:len]
return(as.integer(table(ids)))
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"))
## }

14
app/CMakeLists.txt Normal file
View File

@ -0,0 +1,14 @@
configure_file(poet.h.in poet.h)
if(POET_USE_PRM_BACKEND)
set(poet_SRC poet_prm.cpp)
else()
set(poet_SRC poet.cpp)
endif()
add_executable(poet ${poet_SRC})
target_include_directories(poet PUBLIC "${CMAKE_CURRENT_BINARY_DIR}")
target_link_libraries(poet PUBLIC poet_lib MPI::MPI_CXX)
#target_compile_definitions(poet PRIVATE OMPI_SKIP_MPICXX)
install(TARGETS poet DESTINATION bin)

361
app/poet.cpp Normal file
View File

@ -0,0 +1,361 @@
/*
** 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 "poet/ChemistryModule.hpp"
#include <RInside.h>
#include <Rcpp.h>
#include <Rcpp/internal/wrap.h>
#include <cstdint>
#include <cstdlib>
#include <poet/DiffusionModule.hpp>
#include <poet/Grid.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;
}
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, Grid &grid,
ChemistryParams &chem_params,
const GridParams &g_params, uint32_t nxyz_master) {
DiffusionModule diffusion(poet::DiffusionParams(R), grid);
/* 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,
MPI_COMM_WORLD);
set_chem_parameters(chem, nxyz_master, chem_params.database_path);
chem.RunInitFile(chem_params.input_script);
chem.InitializeField(grid.GetTotalCellCount(), DFToHashMap(g_params.init_df));
if (params.getNumParams().dht_enabled) {
chem.SetDHTEnabled(true, params.getNumParams().dht_size_per_process);
if (!params.getDHTSignifVector().empty()) {
chem.SetDHTSignifVector(params.getDHTSignifVector());
}
if (!params.getDHTPropTypeVector().empty()) {
chem.SetDHTPropTypeVector(params.getDHTPropTypeVector());
}
if (!params.getDHTFile().empty()) {
chem.ReadDHTFile(params.getDHTFile());
}
if (params.getNumParams().dht_snaps > 0) {
chem.SetDHTSnaps(params.getNumParams().dht_snaps, params.getOutDir());
}
}
StateMemory *chem_state = grid.RegisterState("state_C", chem.GetPropNames());
chem_state->mem = chem.GetField();
/* SIMULATION LOOP */
double dStartTime = MPI_Wtime();
for (uint32_t iter = 1; iter < maxiter + 1; iter++) {
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;
/* displaying iteration number, with C++ and R iterator */
cout << "CPP: Going through iteration " << iter << endl;
cout << "CPP: R's $iter: " << ((uint32_t)(R.parseEval("mysetup$iter")))
<< ". Iteration" << endl;
/* run transport */
// TODO: transport to diffusion
diffusion.simulate(dt);
grid.PreModuleFieldCopy(tick++);
cout << "CPP: Chemistry" << endl;
chem.SetTimeStep(dt);
chem.SetField(chem_state->mem);
chem.SetTimeStep(dt);
chem.RunCells();
chem_state->mem = chem.GetField();
grid.WriteFieldsToR(R);
grid.PreModuleFieldCopy(tick++);
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)");
cout << endl
<< "CPP: End of *coupling* iteration " << iter << "/" << maxiter
<< endl
<< endl;
// MPI_Barrier(MPI_COMM_WORLD);
} // 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.getNumParams().dht_enabled) {
R["dht_hits"] = Rcpp::wrap(chem.GetWorkerDHTHits());
R.parseEvalQ("profiling$dht_hits <- dht_hits");
R["dht_miss"] = Rcpp::wrap(chem.GetWorkerDHTMiss());
R.parseEvalQ("profiling$dht_miss <- dht_miss");
R["dht_evictions"] = Rcpp::wrap(chem.GetWorkerDHTEvictions());
R.parseEvalQ("profiling$dht_evictions <- dht_evictions");
R["dht_get_time"] = Rcpp::wrap(chem.GetWorkerDHTGetTimings());
R.parseEvalQ("profiling$dht_get_time <- dht_get_time");
R["dht_fill_time"] = Rcpp::wrap(chem.GetWorkerDHTFillTimings());
R.parseEvalQ("profiling$dht_fill_time <- dht_fill_time");
}
chem.MasterLoopBreak();
diffusion.end();
return MPI_Wtime() - dStartTime;
}
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);
if (world_rank == 0) {
cout << "Running POET in version " << poet_version << endl << endl;
}
if (world_rank > 0) {
{
uint32_t c_size;
MPI_Bcast(&c_size, 1, MPI_UINT32_T, 0, MPI_COMM_WORLD);
char *buffer = new char[c_size + 1];
MPI_Bcast(buffer, c_size, MPI_CHAR, 0, MPI_COMM_WORLD);
buffer[c_size] = '\0';
// ChemistryModule worker(nxyz, nxyz, MPI_COMM_WORLD);
ChemistryModule worker =
poet::ChemistryModule::createWorker(MPI_COMM_WORLD);
set_chem_parameters(worker, worker.GetWPSize(), std::string(buffer));
delete[] buffer;
worker.WorkerLoop();
}
MPI_Barrier(MPI_COMM_WORLD);
cout << "CPP: finished, cleanup of process " << world_rank << endl;
MPI_Finalize();
return EXIT_SUCCESS;
}
/* initialize R runtime */
RInside R(argc, argv);
/*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;
}
cout << "CPP: R Init (RInside) on process " << world_rank << endl;
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);
Grid grid;
GridParams g_params(R);
grid.InitModuleFromParams(g_params);
grid.PushbackModuleFlow(poet::DIFFUSION_MODULE_NAME, CHEMISTRY_MODULE_NAME);
grid.PushbackModuleFlow(CHEMISTRY_MODULE_NAME, poet::DIFFUSION_MODULE_NAME);
params.initVectorParams(R, grid.GetSpeciesCount());
// MDL: store all parameters
if (world_rank == 0) {
cout << "CPP: Calling R Function to store calling parameters" << endl;
R.parseEvalQ("StoreSetup(setup=mysetup)");
}
if (world_rank == 0) {
cout << "CPP: Init done on process with rank " << world_rank << endl;
}
// MPI_Barrier(MPI_COMM_WORLD);
poet::ChemistryParams chem_params(R);
/* THIS IS EXECUTED BY THE MASTER */
std::string db_path = chem_params.database_path;
uint32_t c_size = db_path.size();
MPI_Bcast(&c_size, 1, MPI_UINT32_T, 0, MPI_COMM_WORLD);
MPI_Bcast(db_path.data(), c_size, MPI_CHAR, 0, MPI_COMM_WORLD);
uint32_t nxyz_master = (world_size == 1 ? grid.GetTotalCellCount() : 1);
dSimTime = RunMasterLoop(params, R, grid, chem_params, g_params, nxyz_master);
cout << "CPP: finished simulation loop" << endl;
cout << "CPP: start timing profiling" << endl;
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);
cout << "CPP: Done! Results are stored as R objects into <"
<< params.getOutDir() << "/timings.rds>" << endl;
MPI_Barrier(MPI_COMM_WORLD);
cout << "CPP: finished, cleanup of process " << world_rank << endl;
MPI_Finalize();
if (world_rank == 0) {
cout << "CPP: done, bye!" << endl;
}
exit(0);
}

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

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

287
app/poet_prm.cpp Normal file
View File

@ -0,0 +1,287 @@
/*
** 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 "poet/ChemistryModule.hpp"
#include <RInside.h>
#include <Rcpp.h>
#include <Rcpp/internal/wrap.h>
#include <cstdint>
#include <cstdlib>
#include <poet/DiffusionModule.hpp>
#include <poet/Grid.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;
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, Grid &grid,
ChemistryParams &chem_params,
const GridParams &g_params, uint32_t nxyz_master) {
DiffusionModule diffusion(poet::DiffusionParams(R), grid);
/* 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(grid.GetTotalCellCount(), MPI_COMM_WORLD);
set_chem_parameters(chem, nxyz_master, chem_params.database_path);
chem.RunInitFile(chem_params.input_script);
chem.SetSelectedOutputOn(true);
chem.SetTimeStep(0);
chem.RunCells();
StateMemory *chem_state = grid.RegisterState("state_C", chem.GetPropNames());
auto &chem_field = chem_state->mem;
chem_field = chem.GetField();
/* SIMULATION LOOP */
double dStartTime = MPI_Wtime();
for (uint32_t iter = 1; iter < maxiter + 1; iter++) {
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;
/* displaying iteration number, with C++ and R iterator */
cout << "CPP: Going through iteration " << iter << endl;
cout << "CPP: R's $iter: " << ((uint32_t)(R.parseEval("mysetup$iter")))
<< ". Iteration" << endl;
/* run transport */
// TODO: transport to diffusion
diffusion.simulate(dt);
grid.PreModuleFieldCopy(tick++);
cout << "CPP: Chemistry" << endl;
chem.SetTimeStep(dt);
chem.SetConcentrations(chem_field);
chem.SetTimeStep(dt);
chem.RunCells();
chem_field = chem.GetField();
grid.WriteFieldsToR(R);
grid.PreModuleFieldCopy(tick++);
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)");
cout << endl
<< "CPP: End of *coupling* iteration " << iter << "/" << maxiter
<< endl
<< endl;
// MPI_Barrier(MPI_COMM_WORLD);
} // END SIMULATION LOOP
R.parseEvalQ("profiling <- list()");
R["simtime_chemistry"] = chem.GetChemistryTime();
R.parseEvalQ("profiling$simtime_chemistry <- simtime_chemistry");
chem.MpiWorkerBreak();
diffusion.end();
return MPI_Wtime() - dStartTime;
}
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);
if (world_rank == 0) {
cout << "Running POET in version " << poet_version << endl << endl;
}
if (world_rank > 0) {
{
uint32_t nxyz;
MPI_Bcast(&nxyz, 1, MPI_UINT32_T, 0, MPI_COMM_WORLD);
ChemistryModule worker(nxyz, MPI_COMM_WORLD);
worker.MpiWorker();
}
MPI_Barrier(MPI_COMM_WORLD);
cout << "CPP: finished, cleanup of process " << world_rank << endl;
MPI_Finalize();
return EXIT_SUCCESS;
}
/* initialize R runtime */
RInside R(argc, argv);
/*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;
}
cout << "CPP: R Init (RInside) on process " << world_rank << endl;
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);
Grid grid;
GridParams g_params(R);
grid.InitModuleFromParams(g_params);
grid.PushbackModuleFlow(poet::DIFFUSION_MODULE_NAME, CHEMISTRY_MODULE_NAME);
grid.PushbackModuleFlow(CHEMISTRY_MODULE_NAME, poet::DIFFUSION_MODULE_NAME);
params.initVectorParams(R, grid.GetSpeciesCount());
// MDL: store all parameters
if (world_rank == 0) {
cout << "CPP: Calling R Function to store calling parameters" << endl;
R.parseEvalQ("StoreSetup(setup=mysetup)");
}
if (world_rank == 0) {
cout << "CPP: Init done on process with rank " << world_rank << endl;
}
// MPI_Barrier(MPI_COMM_WORLD);
poet::ChemistryParams chem_params(R);
/* THIS IS EXECUTED BY THE MASTER */
uint32_t nxyz = grid.GetTotalCellCount();
MPI_Bcast(&nxyz, 1, MPI_UINT32_T, 0, MPI_COMM_WORLD);
uint32_t nxyz_master = grid.GetTotalCellCount();
dSimTime = RunMasterLoop(params, R, grid, chem_params, g_params, nxyz_master);
cout << "CPP: finished simulation loop" << endl;
cout << "CPP: start timing profiling" << endl;
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);
cout << "CPP: Done! Results are stored as R objects into <"
<< params.getOutDir() << "/timings.rds>" << endl;
MPI_Barrier(MPI_COMM_WORLD);
cout << "CPP: finished, cleanup of process " << world_rank << endl;
MPI_Finalize();
if (world_rank == 0) {
cout << "CPP: done, bye!" << endl;
}
exit(0);
}

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

20
bench/barite/CMakeLists.txt
View File

@ -1,20 +0,0 @@
# Create a list of files
set(bench_files
barite_200.R
barite_het.R
)
set(runtime_files
barite_200_rt.R
barite_het_rt.R
)
# add_custom_target(barite_bench DEPENDS ${bench_files} ${runtime_files})
ADD_BENCH_TARGET(barite_bench
bench_files
runtime_files
"barite"
)
add_dependencies(${BENCHTARGET} barite_bench)

134
bench/barite/README.org
View File

@ -1,134 +0,0 @@
#+TITLE: Description of =barite= benchmark
#+AUTHOR: MDL <delucia@gfz-potsdam.de>
#+DATE: 2023-08-26
#+STARTUP: inlineimages
#+LATEX_CLASS_OPTIONS: [a4paper,9pt]
#+LATEX_HEADER: \usepackage{fullpage}
#+LATEX_HEADER: \usepackage{amsmath, systeme}
#+LATEX_HEADER: \usepackage{graphicx}
#+LATEX_HEADER: \usepackage{charter}
#+OPTIONS: toc:nil
* Quick start
#+begin_src sh :language sh :frame single
mpirun -np 4 ./poet barite.R barite_results
mpirun -np 4 ./poet --interp barite_interp_eval.R barite_results
#+end_src
* List of Files
- =barite_het.R=: POET input script with homogeneous zones for a 5x2 simulation
grid
- =barite_200.R=: POET input script for a 200x200 simulation
grid
- =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
* Chemical system
The benchmark depicts an isotherm porous system at *25 °C* where pure
water is initially at equilibrium with *celestite* (strontium sulfate;
brute formula: SrSO_{4}). A solution containing only dissolved Ba^{2+}
and Cl^- diffuses into the system causing celestite dissolution. The
increased concentration of dissolved sulfate SO_{4}^{2-} induces
precipitation of *barite* (barium sulfate; brute formula:
BaSO_{4}^{2-}). The overall reaction can be written:
Ba^{2+} + celestite \rightarrow barite + Sr^{2+}
Both celestite dissolution and barite precipitation are calculated
using a kinetics rate law based on transition state theory:
rate = -S_{m} k_{barite} (1-SR_{m})
where the reaction rate has units mol/s, S_{m} (m^{2}) is the reactive
surface area, k (mol/m^{2}/s) is the kinetic coefficient, and SR is
the saturation ratio, i.e., the ratio of the ion activity product of
the reacting species and the solubility constant, calculated
internally by PHREEQC from the speciated solution.
For barite, the reaction rate is computed as sum of two mechanisms,
r_{/acid/} and r_{/neutral/}:
rate_{barite} = S_{barite} (k_{/acid/} + k_{/neutral/}) * (1 - SR_{barite})
where:
k_{/acid/} = 10^{-6.9} e^{-30800 / R} \cdot act(H^{+})^{0.22}
k_{/neutral/} = 10^{-7.9} e^{-30800 / R}
R (8.314462 J K^{-1} mol^{-1}) is the gas constant.
For celestite the kinetic law considers only the acidic mechanism and
reads:
rate_{celestite} = S_{celestite} 10^{-5.66} e^{-23800 / R} \cdot
act(H^{+})^{0.109} \cdot (1 - SR_{celestite})
The kinetic rates as implemented in the =db_barite.dat= file accepts
one parameter which represents reactive surface area in m^{2}. For the
benchmarks the surface areas are set to
- S_{barite}: 50 m^{2}
- S_{celestite}: 10 m^{2}
A starting seed for barite is given at 0.001 mol in each domain
element.
* Nucleation (TODO)
Geochemical benchmark inspired by Tranter et al. (2021) without
nucleation.
* POET simulations
Currently these benchmarks are pure diffusion simulations. There are 7
transported species: H, O, Charge, Ba, Cl, S(6), Sr.
** =barite.R=
- Grid discretization: square domain of 1 \cdot 1 m^{2} discretized in
20x20 cells
- Boundary conditions: E, S and W sides of the domain are closed; the
N boundary has a *fixed concentration* (Dirichlet) of 0.1 molal
BaCl_{2}
- Diffusion coefficients: isotropic homogeneous \alpha = 1E-06
- Time steps & iterations: 20 iteration with \Delta t = 250 s
- *DHT* parameters:
| H | O | Charge | Ba | Cl | S(6) | Sr |
| 10 | 10 | 3 | 5 | 5 | 5 | 5 |
** =barite_interp_eval.R=
- Grid discretization: rectangular domain of 40 \cdot 20 m^{2}
discretized in 400 \cdot 200 cells
- Boundary conditions: all boundaries are closed. The center of the
domain at indeces (200, 100) has fixed concentration of 0.1 molal of
BaCl_{2}
- Diffusion coefficients: isotropic homogeneous \alpha = 1E-06
- Time steps & iterations: 200 iterations with \Delta t = 250 s
- *DHT* parameters:
| H | O | Charge | Ba | Cl | S(6) | Sr |
| 10 | 10 | 3 | 5 | 5 | 5 | 5 |
* References
- Tranter, Wetzel, De Lucia and Kühn (2021): Reactive transport model
of kinetically controlled celestite to barite replacement, Advances
in Geosciences, 56, 57-–65, 2021.
https://doi.org/10.5194/adgeo-56-57-20211

32
bench/barite/barite.pqi
View File

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

195
bench/barite/db_barite.dat
View File

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

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.

18
bench/dolo/CMakeLists.txt
View File

@ -1,18 +0,0 @@
set(bench_files
dolo_inner_large.R
dolo_interp.R
)
set(runtime_files
dolo_inner_large_rt.R
dolo_interp_rt.R
)
ADD_BENCH_TARGET(
dolo_bench
bench_files
runtime_files
"dolo"
)
add_dependencies(${BENCHTARGET} dolo_bench)

159
bench/dolo/README.org
View File

@ -1,159 +0,0 @@
#+TITLE: Description of =dolo= benchmark
#+AUTHOR: MDL <delucia@gfz-potsdam.de>
#+DATE: 2023-08-26
#+STARTUP: inlineimages
#+LATEX_CLASS_OPTIONS: [a4paper,9pt]
#+LATEX_HEADER: \usepackage{fullpage}
#+LATEX_HEADER: \usepackage{amsmath, systeme}
#+LATEX_HEADER: \usepackage{graphicx}
#+LATEX_HEADER: \usepackage{charter}
#+OPTIONS: toc:nil
* Quick start
#+begin_src sh :language sh :frame single
mpirun -np 4 ./poet dolo_diffu_inner.R dolo_diffu_inner_res
mpirun -np 4 ./poet --dht --interp dolo_interp_long.R dolo_interp_long_res
#+end_src
* List of Files
- =dolo_interp.R=: POET input script for a 400x200 simulation
grid
- =dolo_diffu_inner_large.R=: POET input script for a 400x200
simulation grid
- =phreeqc_kin.dat=: PHREEQC database containing the kinetic expressions
for dolomite and celestite, stripped down from =phreeqc.dat=
- =dol.pqi=: PHREEQC input script for the chemical system
* Chemical system
This model describes a simplified version of /dolomitization/ of
calcite, itself a complex and not yet fully understood natural process
which is observed naturally at higher temperatures (Möller and De
Lucia, 2020). Variants of such model have been widely used in many
instances especially for the purpose of benchmarking numerical codes
(De Lucia et al., 2021 and references therein).
We consider an isotherm porous system at *25 °C* in which pure water
is initially at equilibrium with *calcite* (calcium carbonate; brute
formula: CaCO_{3}). An MgCl_{2} solution enters the system causing
calcite dissolution. The released excess concentration of dissolved
calcium Ca^{2+} and carbonate (CO_{3}^{2-}) induces after a while
supersaturation and hence precipitation of *dolomite*
(calcium-magnesium carbonate; brute formula: CaMg(CO_{3})_{2}). The
overall /dolomitization/ reaction can be written:
Mg^{2+} + 2 \cdot calcite \rightarrow dolomite + 2 \cdot Ca^{2+}
The precipitation of dolomite continues until enough Ca^{2+} is
present in solution. Further injection of MgCl_{2} changes its
saturation state causing its dissolution too. After enough time, the
whole system has depleted all minerals and the injected MgCl_{2}
solution fills up the domain.
Both calcite dissolution and dolomite precipitation/dissolution follow
a kinetics rate law based on transition state theory (Palandri and
Karhaka, 2004; De Lucia et al., 2021).
rate = -S_{m} k_{m} (1-SR_{m})
where the reaction rate has units mol/s, S_{m} (m^{2}) is the reactive
surface area, k_{m} (mol/m^{2}/s) is the kinetic coefficient, and SR
is the saturation ratio, i.e., the ratio of the ion activity product
of the reacting species and the solubility constant, calculated
internally by PHREEQC from the speciated solution.
For dolomite, the kinetic coefficient results from the sum of two
mechanisms, r_{/acid/} and r_{/neutral/}:
rate_{dolomite} = S_{dolomite} (k_{/acid/} + k_{/neutral/}) * (1 - SR_{dolomite})
where:
k_{/acid/} = 10^{-3.19} e^{-36100 / R} \cdot act(H^{+})^{0.5}
k_{/neutral/} = 10^{-7.53} e^{-52200 / R}
R (8.314462 J K^{-1} mol^{-1}) is the gas constant.
Similarly, the kinetic law for calcite reads:
k_{/acid/} = 10^{-0.3} e^{-14400 / R} \cdot act(H^{+})^{0.5}
k_{/neutral/} = 10^{-5.81} e^{-23500 / R}
The kinetic laws as implemented in the =phreeqc_kin.dat= file accepts
one parameter which represents reactive surface area in m^{2}. For the
benchmarks the surface areas are set to
- S_{dolomite}: 0.005 m^{2}
- S_{calcite}: 0.05 m^{2}
The initial content of calcite in the domain is of 0.0002 mol per kg
of water. A constant partial pressure of 10^{-1675} atm of O_{2(g)} is
maintained at any time in the domain in order to fix the redox
potential of the solution to an oxidizing state (pe around 9).
Note that Cl is unreactive in this system and only effects the
computation of the activities in solution.
* POET simulations
Several benchmarks based on the same chemical system are defined here
with different grid sizes, resolution and boundary conditions. The
transported elemental concentrations are 7: C(4), Ca, Cl, Mg and the
implicit total H, total O and Charge as required by PHREEQC_RM.
** =dolo_diffu_inner.R=
- Grid discretization: square domain of 1 \cdot 1 m^{2} discretized in
100x100 cells
- Boundary conditions: All sides of the domain are closed. *Fixed
concentration* of 0.001 molal of MgCl_{2} is defined in the domain
cell (20, 20) and of 0.002 molal MgCl_{2} at cells (60, 60) and
(80, 80)
- Diffusion coefficients: isotropic homogeneous \alpha = 1E-06
- Time steps & iterations: 10 iterations with \Delta t of 200 s
- *DHT* parameters:
| H | O | Charge | C(4) | Ca | Cl | Mg | Calcite | Dolomite |
| 10 | 10 | 3 | 5 | 5 | 5 | 5 | 5 | 5 |
- Hooks: the following hooks are defined:
1. =dht_fill=:
2. =dht_fuzz=:
3. =interp_pre_func=:
4. =interp_post_func=:
** =dolo_interp_long.R=
- Grid discretization: rectangular domain of 5 \cdot 2.5 m^{2}
discretized in 400 \times 200 cells
- Boundary conditions: *Fixed concentrations* equal to the initial
state are imposed at all four sides of the domain. *Fixed
concentration* of 0.001 molal of MgCl_{2} is defined in the domain
center at cell (100, 50)
- Diffusion coefficients: isotropic homogeneous \alpha = 1E-06
- Time steps & iterations: 20000 iterations with \Delta t of 200 s
- *DHT* parameters:
| H | O | Charge | C(4) | Ca | Cl | Mg | Calcite | Dolomite |
| 10 | 10 | 3 | 5 | 5 | 5 | 5 | 5 | 5 |
- Hooks: the following hooks are defined:
1. =dht_fill=:
2. =dht_fuzz=:
3. =interp_pre_func=:
4. =interp_post_func=:
* References
- De Lucia, Kühn, Lindemann, Lübke, Schnor: POET (v0.1): speedup of
many-core parallel reactive transport simulations with fast DHT
lookups, Geosci. Model Dev., 14, 7391–7409, 2021.
https://doi.org/10.5194/gmd-14-7391-2021
- Möller, Marco De Lucia: The impact of Mg^{2+} ions on equilibration
of Mg-Ca carbonates in groundwater and brines, Geochemistry
80, 2020. https://doi.org/10.1016/j.chemer.2020.125611
- Palandri, Kharaka: A Compilation of Rate Parameters of Water-Mineral
Interaction Kinetics for Application to Geochemical Modeling, Report
2004-1068, USGS, 2004.

43
bench/dolo/dol.pqi
View File

@ -1,43 +0,0 @@
SOLUTION 1
units mol/kgw
water 1
temperature 25
pH 7
pe 4
PURE 1
Calcite 0.0 1
END
RUN_CELLS
-cells 1
COPY solution 1 2
PURE 2
O2g -0.1675 10
KINETICS 2
Calcite
-m 0.000207
-parms 0.05
-tol 1e-10
Dolomite
-m 0.0
-parms 0.005
-tol 1e-10
END
SOLUTION 3
units mol/kgw
water 1
temp 25
Mg 0.001
Cl 0.002
END
SOLUTION 4
units mol/kgw
water 1
temp 25
Mg 0.002
Cl 0.004
END

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
)

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
)

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

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

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

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

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

@ -0,0 +1,146 @@
## Time-stamp: "Last modified 2023-01-10 13:51:40 delucia"
database <- normalizePath("../share/poet/examples/phreeqc_kin.dat")
input_script <- normalizePath("../share/poet/bench/dolo_inner.pqi")
#################################################################
## Section 1 ##
## Grid initialization ##
#################################################################
n <- 100
m <- 100
types <- c("scratch", "phreeqc", "rds")
init_cell <- list(
"H" = 110.683,
"O" = 55.3413,
"Charge" = -5.0822e-19,
"C" = 1.2279E-4,
"Ca" = 1.2279E-4,
"Cl" = 0,
"Mg" = 0,
"O2g" = 0.499957,
"Calcite" = 2.07e-4,
"Dolomite" = 0
)
grid <- list(
n_cells = c(n, m),
s_cells = c(1, 1),
type = types[1],
init_cell = as.data.frame(init_cell),
props = names(init_cell),
database = database,
input_script = input_script
)
##################################################################
## Section 2 ##
## Diffusion parameters and boundary conditions ##
##################################################################
## initial conditions
init_diffu <- c(
"H" = 110.683,
"O" = 55.3413,
"Charge" = -5.0822e-19,
"C" = 1.2279E-4,
"Ca" = 1.2279E-4,
"Cl" = 0,
"Mg" = 0
)
## diffusion coefficients
alpha_diffu <- c(
"H" = 1E-6,
"O" = 1E-6,
"Charge" = 1E-6,
"C" = 1E-6,
"Ca" = 1E-6,
"Cl" = 1E-6,
"Mg" = 1E-6
)
## list of boundary conditions/inner nodes
vecinj_diffu <- list(
list(
"H" = 110.683,
"O" = 55.3413,
"Charge" = 1.90431e-16,
"C" = 0,
"Ca" = 0,
"Cl" = 0.002,
"Mg" = 0.001
),
list(
"H" = 110.683,
"O" = 55.3413,
"Charge" = 1.90431e-16,
"C" = 0,
"Ca" = 0.0,
"Cl" = 0.004,
"Mg" = 0.002
)
)
vecinj_inner <- list(
l1 = c(1,20,20),
l2 = c(2,80,80),
l3 = c(2,60,80)
)
boundary <- list(
# "N" = c(1, rep(0, n-1)),
"N" = rep(0, n),
"E" = rep(0, n),
"S" = rep(0, n),
"W" = rep(0, n)
)
diffu_list <- names(alpha_diffu)
diffusion <- list(
init = init_diffu,
vecinj = do.call(rbind.data.frame, vecinj_diffu),
vecinj_inner = vecinj_inner,
vecinj_index = boundary,
alpha = alpha_diffu
)
#################################################################
## Section 3 ##
## Chemistry module (Phreeqc) ##
#################################################################
## # Needed when using DHT
signif_vector <- c(10, 10, 2, 5, 5, 5, 5, 0, 5, 5)
prop_type <- c("", "", "", "act", "act", "act", "act", "ignore", "", "")
prop <- names(init_cell)
chemistry <- list(
database = database,
input_script = input_script
)
#################################################################
## Section 4 ##
## Putting all those things together ##
#################################################################
iterations <- 1000
dt <- 200
setup <- list(
grid = grid,
diffusion = diffusion,
chemistry = chemistry,
iterations = iterations,
timesteps = rep(dt, iterations),
store_result = TRUE,
out_save = c(5, iterations)
)

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

@ -0,0 +1,146 @@
## Time-stamp: "Last modified 2023-01-10 13:51:40 delucia"
database <- normalizePath("../share/poet/examples/phreeqc_kin.dat")
input_script <- normalizePath("../share/poet/bench/dolo_inner.pqi")
#################################################################
## Section 1 ##
## Grid initialization ##
#################################################################
n <- 2000
m <- 1000
types <- c("scratch", "phreeqc", "rds")
init_cell <- list(
"H" = 110.683,
"O" = 55.3413,
"Charge" = -5.0822e-19,
"C" = 1.2279E-4,
"Ca" = 1.2279E-4,
"Cl" = 0,
"Mg" = 0,
"O2g" = 0.499957,
"Calcite" = 2.07e-4,
"Dolomite" = 0
)
grid <- list(
n_cells = c(n, m),
s_cells = c(2, 1),
type = types[1],
init_cell = as.data.frame(init_cell),
props = names(init_cell),
database = database,
input_script = input_script
)
##################################################################
## Section 2 ##
## Diffusion parameters and boundary conditions ##
##################################################################
## initial conditions
init_diffu <- c(
"H" = 110.683,
"O" = 55.3413,
"Charge" = -5.0822e-19,
"C" = 1.2279E-4,
"Ca" = 1.2279E-4,
"Cl" = 0,
"Mg" = 0
)
## diffusion coefficients
alpha_diffu <- c(
"H" = 1E-6,
"O" = 1E-6,
"Charge" = 1E-6,
"C" = 1E-6,
"Ca" = 1E-6,
"Cl" = 1E-6,
"Mg" = 1E-6
)
## list of boundary conditions/inner nodes
vecinj_diffu <- list(
list(
"H" = 110.683,
"O" = 55.3413,
"Charge" = 1.90431e-16,
"C" = 0,
"Ca" = 0,
"Cl" = 0.002,
"Mg" = 0.001
),
list(
"H" = 110.683,
"O" = 55.3413,
"Charge" = 1.90431e-16,
"C" = 0,
"Ca" = 0.0,
"Cl" = 0.004,
"Mg" = 0.002
)
)
vecinj_inner <- list(
l1 = c(1,400,200),
l2 = c(2,1400,800),
l3 = c(2,1600,800)
)
boundary <- list(
# "N" = c(1, rep(0, n-1)),
"N" = rep(0, n),
"E" = rep(0, m),
"S" = rep(0, n),
"W" = rep(0, m)
)
diffu_list <- names(alpha_diffu)
diffusion <- list(
init = init_diffu,
vecinj = do.call(rbind.data.frame, vecinj_diffu),
vecinj_inner = vecinj_inner,
vecinj_index = boundary,
alpha = alpha_diffu
)
#################################################################
## Section 3 ##
## Chemistry module (Phreeqc) ##
#################################################################
## # Needed when using DHT
signif_vector <- c(10, 10, 2, 5, 5, 5, 5, 0, 5, 5)
prop_type <- c("", "", "", "act", "act", "act", "act", "ignore", "", "")
prop <- names(init_cell)
chemistry <- list(
database = database,
input_script = input_script
)
#################################################################
## Section 4 ##
## Putting all those things together ##
#################################################################
iterations <- 500
dt <- 50
setup <- list(
grid = grid,
diffusion = diffusion,
chemistry = chemistry,
iterations = iterations,
timesteps = rep(dt, iterations),
store_result = TRUE,
out_save = c(5, iterations)
)

28
bench/dolo_diffu_inner/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/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
)

20
bench/surfex/CMakeLists.txt
View File

@ -1,20 +0,0 @@
set(bench_files
# surfex.R
# ex.R
PoetEGU_surfex_500.R
)
set(runtime_files
# surfex_rt.R
# ex_rt.R
PoetEGU_surfex_500_rt.R
)
ADD_BENCH_TARGET(
surfex_bench
bench_files
runtime_files
"surfex"
)
add_dependencies(${BENCHTARGET} surfex_bench)

63
bench/surfex/ExBase.pqi
View File

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

100
bench/surfex/README.org
View File

@ -1,100 +0,0 @@
#+TITLE: Description of =surfex= benchmark
#+AUTHOR: MDL <delucia@gfz-potsdam.de>
#+DATE: 2023-08-26
#+STARTUP: inlineimages
#+LATEX_CLASS_OPTIONS: [a4paper,9pt]
#+LATEX_HEADER: \usepackage{fullpage}
#+LATEX_HEADER: \usepackage{amsmath, systeme}
#+LATEX_HEADER: \usepackage{graphicx}
#+LATEX_HEADER: \usepackage{charter}
#+OPTIONS: toc:nil
* Quick start
#+begin_src sh :language sh :frame single
mpirun -np 4 ./poet ex.R ex_res
mpirun -np 4 ./poet surfex.R surfex_res
#+end_src
* List of Files
- =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
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
parametrized data for Surface and Exchange, based on the SMILE
Thermodynamic Database (Version 01-November-2021)
* Chemical system
This model describes migration of Uranium radionuclide in Opalinus
clay subject to surface complexation and cation exchange on the
surface of clay minerals. These two processes account for the binding
of aqueous complexes to the surfaces of minerals, which may have a
significant impact on safety of underground nuclear waste repository.
Namely, they can act as retardation buffer for uranium complexes
entering into a natural system. The system is kindly provided by Dr.
T. Hennig and is inspired to the sandy facies BWS-A3 sample from the
Mont Terri underground lab (Hennig and Kühn, 2021).
This chemical system is highly redox-sensitive, and several elements
are defined in significant amounts in different valence states. In
total, 20 elemental concentrations and valences are transported:
C(-4), C(4), Ca, Cl, Fe(2), Fe(3), K, Mg, Na, S(-2), S(2), S(4), S(6),
Sr , U(4), U(5), U(6); plus the total H, total O and Charge implicitly
required by PHREEQC_RM.
** Exchange
The SMILE database defines thermodynamical data for exchange of all
major cations and uranyl-ions on Illite and Montmorillonite. In
PHREEQC terms:
- *Y* for Montmorillonite, with a total amount of 1.2585
milliequivalents and
- *Z* for Illite, with a total amount of 0.9418 meq
** Surface
Here we consider a Donnan diffuse double layer of 0.49 nm. Six
distinct sorption sites are defined:
- Kln_aOH (aluminol site) and Kln_siOH (silanol) for Kaolinite
- For Illite, strong and weak sites Ill_sOH and Ill_wOH respectively
- For Montmorillonite, strong and weak sites Mll_sOH and Mll_wOH
respectively
Refer to the =SurfExBase.pqi= script for the actual numerical values
of the parameters.
* POET simulations
** =ex.R=
This benchmark only considers EXCHANGE, no mineral or SURFACE
complexation is involved.
- Grid discretization: square domain of 1 \cdot 1 m^{2} discretized in
100x100 cells
- Boundary conditions: E, S and W sides of the domain are closed.
*Fixed concentrations* are fixed at the N boundary.
- Diffusion coefficients: isotropic homogeneous \alpha = 1E-06
- Time steps & iterations: 10 iterations with \Delta t of 200 s
- *DHT* is not implemented as of yet for models including SURFACE and
EXCHANGE geochemical processes *TODO*
- Hooks: no hooks defined *TODO*
** =surfex.R=
- Grid discretization: rectangular domain of 1 \cdot 1 m^{2}
discretized in 10 \times 10 cells
- Boundary conditions: E, S and W sides of the domain are closed.
*Fixed concentrations* are fixed at the N boundary.
- Diffusion coefficients: isotropic homogeneous \alpha = 1E-06
- Time steps & iterations: 10 iterations with \Delta t of 200 s
* References
- Hennig, T.; Kühn, M.Surrogate Model for Multi-Component Diffusion of
Uranium through Opalinus Clay on the Host Rock Scale. Appl. Sci.
2021, 11, 786. https://doi.org/10.3390/app11020786

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

37
bench/surfex/ex.R
View File

@ -1,37 +0,0 @@
rows <- 100
cols <- 100
grid_def <- matrix(1, nrow = rows, ncol = cols)
# Define grid configuration for POET model
grid_setup <- list(
pqc_in_file = "./SurfExBase.pqi",
pqc_db_file = "./SMILE_2021_11_01_TH.dat", # Path to the database file for Phreeqc
grid_def = grid_def, # Definition of the grid, containing IDs according to the Phreeqc input script
grid_size = c(1, 1), # Size of the grid in meters
constant_cells = c() # IDs of cells with constant concentration
)
bound_def <- list(
"type" = rep("constant", cols),
"sol_id" = rep(2, cols),
"cell" = seq(1, cols)
)
diffusion_setup <- list(
boundaries = list(
"N" = bound_def
),
alpha_x = 1e-6,
alpha_y = 1e-6
)
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
)

7
bench/surfex/ex_rt.R
View File

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

161
bench/surfex/surfex.R
View File

@ -1,37 +1,142 @@
rows <- 1000
cols <- 1000
## Time-stamp: "Last modified 2023-02-27 18:33:30 delucia"
grid_def <- matrix(1, nrow = rows, ncol = cols)
database <- normalizePath("./SMILE_2021_11_01_TH.dat")
input_script <- normalizePath("./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",
init_cell = as.data.frame(init_cell, check.names = FALSE),
props = names(init_cell),
database = database,
input_script = input_script
)
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)
diffusion <- list(
init = init_cell,
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
data/CMakeLists.txt Normal file
View File

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

180
data/SimDol1D_diffu.R Normal file
View File

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

213
data/SimDol2D_diffu.R Normal file
View File

@ -0,0 +1,213 @@
database <- normalizePath("../share/poet/examples/phreeqc_kin.dat")
input_script <- normalizePath("../share/poet/examples/dol.pqi")
#################################################################
## Section 1 ##
## Grid initialization ##
#################################################################
n <- 100
m <- 100
types <- c("scratch", "phreeqc", "rds")
# initsim <- c("SELECTED_OUTPUT", "-high_precision true",
# "-reset false",
# "USER_PUNCH",
# "-head C Ca Cl Mg pH pe O2g Calcite Dolomite",
# "10 PUNCH TOT(\"C\"), TOT(\"Ca\"), TOT(\"Cl\"), TOT(\"Mg\"), -LA(\"H+\"), -LA(\"e-\"), EQUI(\"O2g\"), EQUI(\"Calcite\"), EQUI(\"Dolomite\")",
# "SOLUTION 1",
# "units mol/kgw",
# "temp 25.0", "water 1",
# "pH 9.91 charge",
# "pe 4.0",
# "C 1.2279E-04",
# "Ca 1.2279E-04",
# "Cl 1E-12",
# "Mg 1E-12",
# "PURE 1",
# "O2g -0.6788 10.0",
# "Calcite 0.0 2.07E-3",
# "Dolomite 0.0 0.0",
# "END")
# needed if init type is set to "scratch"
# prop <- c("C", "Ca", "Cl", "Mg", "pH", "pe", "O2g", "Calcite", "Dolomite")
init_cell <- list(
"H" = 110.683,
"O" = 55.3413,
"Charge" = -5.0822e-19,
"C" = 1.2279E-4,
"Ca" = 1.2279E-4,
"Cl" = 0,
"Mg" = 0,
"O2g" = 0.499957,
"Calcite" = 2.07e-4,
"Dolomite" = 0
)
grid <- list(
n_cells = c(n, m),
s_cells = c(n, m),
type = types[1],
init_cell = as.data.frame(init_cell),
props = names(init_cell),
database = database,
input_script = input_script
)
##################################################################
## Section 2 ##
## Diffusion parameters and boundary conditions ##
##################################################################
init_diffu <- c(
"H" = 110.683,
"O" = 55.3413,
"Charge" = -5.0822e-19,
"C" = 1.2279E-4,
"Ca" = 1.2279E-4,
"Cl" = 0,
"Mg" = 0
)
alpha_diffu <- c(
"H" = 1E-4,
"O" = 1E-4,
"Charge" = 1E-4,
"C" = 1E-4,
"Ca" = 1E-4,
"Cl" = 1E-4,
"Mg" = 1E-4
)
vecinj_diffu <- list(
list(
"H" = 110.683,
"O" = 55.3413,
"Charge" = 1.90431e-16,
"C" = 0,
"Ca" = 0,
"Cl" = 0.002,
"Mg" = 0.001
)
)
#inner_index <- c(5, 15, 25)
#inner_vecinj_index <- rep(1, 3)
#
#vecinj_inner <- cbind(inner_index, inner_vecinj_index)
vecinj_inner <- list(
l1 = c(1,2,2)
)
boundary <- list(
"N" = rep(0, n),
"E" = rep(0, n),
"S" = rep(0, n),
"W" = rep(0, n)
)
diffu_list <- names(alpha_diffu)
diffusion <- list(
init = init_diffu,
vecinj = do.call(rbind.data.frame, vecinj_diffu),
vecinj_inner = vecinj_inner,
vecinj_index = boundary,
alpha = alpha_diffu
)
#################################################################
## Section 3 ##
## Chemitry module (Phreeqc) ##
#################################################################
# db <- RPhreeFile(system.file("extdata", "phreeqc_kin.dat",
# package = "RedModRphree"
# ), is.db = TRUE)
#
# phreeqc::phrLoadDatabaseString(db)
# NOTE: This won't be needed in the future either. Could also be done in a. pqi
# file
base <- c(
"SOLUTION 1",
"units mol/kgw",
"temp 25.0",
"water 1",
"pH 9.91 charge",
"pe 4.0",
"C 1.2279E-04",
"Ca 1.2279E-04",
"Mg 0.001",
"Cl 0.002",
"PURE 1",
"O2g -0.1675 10",
"KINETICS 1",
"-steps 100",
"-step_divide 100",
"-bad_step_max 2000",
"Calcite", "-m 0.000207",
"-parms 0.0032",
"Dolomite",
"-m 0.0",
"-parms 0.00032",
"END"
)
selout <- c(
"SELECTED_OUTPUT", "-high_precision true", "-reset false",
"-time", "-soln", "-temperature true", "-water true",
"-pH", "-pe", "-totals C Ca Cl Mg",
"-kinetic_reactants Calcite Dolomite", "-equilibrium O2g"
)
# Needed when using DHT
signif_vector <- c(10, 10, 10, 7, 7, 7, 7, 0, 5, 5)
prop_type <- c("", "", "", "act", "act", "act", "act", "ignore", "", "")
prop <- names(init_cell)
chemistry <- list(
database = database,
input_script = input_script
)
#################################################################
## Section 4 ##
## Putting all those things together ##
#################################################################
iterations <- 10
setup <- list(
# bound = myboundmat,
base = base,
first = selout,
# initsim = initsim,
# Cf = 1,
grid = grid,
diffusion = diffusion,
chemistry = chemistry,
prop = prop,
immobile = c(7, 8, 9),
kin = c(8, 9),
ann = list(O2g = -0.1675),
# phase = "FLUX1",
# density = "DEN1",
reduce = FALSE,
# snapshots = demodir, ## directory where we will read MUFITS SUM files
# gridfile = paste0(demodir, "/d2ascii.run.GRID.SUM")
# init = init,
# vecinj = vecinj,
# cinj = c(0,1),
# boundary = boundary,
# injections = FALSE,
iterations = iterations,
timesteps = rep(10, iterations)
)

13
data/SimDol2D_diffu.yaml Normal file
View File

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

35
data/dol.pqi Normal file
View File

@ -0,0 +1,35 @@
SELECTED_OUTPUT
-high_precision true
-reset false
-time
-soln
-temperature true
-water true
-pH
-pe
-totals C Ca Cl Mg
-kinetic_reactants Calcite Dolomite
-equilibrium O2g
SOLUTION 1
units mol/kgw
temp 25.0
water 1
pH 9.91 charge
pe 4.0
C 1.2279E-04
Ca 1.2279E-04
Cl 1E-12
Mg 1E-12
PURE 1
O2g -0.1675 10
KINETICS 1
Calcite
-m 0.000207
-parms 0.0032
-tol 1e-10
Dolomite
-m 0.0
-parms 0.00032
-tol 1e-10
END

6
bench/dolo/phreeqc_kin.dat → data/phreeqc_kin.dat
View File

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

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

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()

95
docs/Input_Scripts.md Normal file
View File

@ -0,0 +1,95 @@
# 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* |
| `init_cell` | Data Frame | Containing all exactly one value per species to initialize the field. |
| `props` | String Vector | Names of all species |
| `database` | String | Path to Phreeqc database |
| `input_script` | String | Path to Phreeqc initial script |
## Diffusion parameters
| 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 |
## 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 |
### DHT setup
| name | type | description |
|-----------------|----------------|---------------------------------------------------------------------------------|
| `signif_vector` | Numeric Vector | Indicates significant digits to use for DHT rounding. Order of `props` vector. |
| `prop_type` | String Vector | Set type of species for rounding, can be left blank or set to *act* or *ignore* |

60
docs/Output.md
View File

@ -35,50 +35,34 @@ corresponding values can be found in `<OUTPUT_DIRECTORY>/timings.rds`
and possible to read out within a R runtime with
`readRDS("timings.rds")`. There you will find the following values:
| Value | Description |
| --------- | -------------------------------------------------------------------------- |
| simtime | time spent in whole simulation loop without any initialization and cleanup |
| chemistry | measured time in *chemistry* subroutine |
| diffusion | measured time in *diffusion* subroutine |
| Value | Description |
|--------------------|----------------------------------------------------------------------------|
| simtime | time spent in whole simulation loop without any initialization and cleanup |
| simtime\_transport | measured time in *transport* subroutine |
| simtime\_chemistry | measured time in *chemistry* subroutine (actual parallelized part) |
### Chemistry subsetting
### chemistry subsetting
| Value | Description |
| ------------- | --------------------------------------------------------- |
| simtime | overall runtime of chemistry |
| loop | time spent in send/recv loop of master |
| sequential | sequential part of the master (e.g. shuffling field) |
| idle\_master | idling time of the master waiting for workers |
| idle\_worker | idling time (waiting for work from master) of the workers |
| phreeqc\_time | accumulated times for Phreeqc calls of every worker |
If running parallel there are also measured timings which are subsets of
*simtime\_​chemistry*.
#### DHT usage
| Value | Description |
|-----------------------|-----------------------------------------------------------|
| chemistry\_loop | time spent in send/​recv loop of master |
| chemistry\_sequential | sequential part of master chemistry |
| idle\_master | idling time (waiting for any free worker) of the master |
| idle\_worker | idling time (waiting for work from master) of the workers |
| phreeqc\_time | accumulated times for Phreeqc calls of every worker |
### DHT usage {#DHT-usage}
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 8fdfd113dcb4ad1a31705ff8ddb13a21a505bad8

1
ext/litephreeqc

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

1
ext/phreeqcrm Submodule

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

2
ext/tug

@ -1 +1 @@
Subproject commit 9c4aeee410c71d064f7567143d4f8d6451ade75a
Subproject commit 79d7a32fc25fd97b78f320bba854b1db64e059f7

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

@ -1,24 +1,15 @@
#ifndef CHEMISTRYMODULE_H_
#define CHEMISTRYMODULE_H_
#include "DataStructures/Field.hpp"
#include "DataStructures/NamedVector.hpp"
#include "IrmResult.h"
#include "PhreeqcRM.h"
#include "poet/DHT_Wrapper.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>
#include <mpi.h>
#include <string>
#include <unordered_map>
#include <vector>
namespace poet {
@ -26,8 +17,22 @@ namespace poet {
* \brief Wrapper around PhreeqcRM to provide POET specific parallelization with
* easy access.
*/
class ChemistryModule {
class ChemistryModule : public PhreeqcRM {
public:
#ifdef POET_USE_PRM
/**
* Creates a new instance of Chemistry module with given grid cell count and
* using MPI communicator.
*
* Constructor is just a wrapper around PhreeqcRM's constructor, so just calls
* the base class constructor.
*
* \param nxyz Count of grid cells.
* \param communicator MPI communicator where work will be distributed.
*/
ChemistryModule(uint32_t nxyz, MPI_Comm communicator)
: PhreeqcRM(nxyz, communicator), n_cells(nxyz){};
#else
/**
* Creates a new instance of Chemistry module with given grid cell count, work
* package size and communicator.
@ -43,80 +48,62 @@ 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, MPI_Comm communicator);
/**
* Deconstructor, which frees DHT data structure if used.
*/
~ChemistryModule();
#endif
/**
* Parses the input script and extract information needed during runtime.
*
* **Only run by master**.
*
* Database must be loaded beforehand.
*
* \param input_script_path Path to input script to parse.
*/
void RunInitFile(const std::string &input_script_path);
void masterSetField(Field field);
/**
* Run the chemical simulation with parameters set.
*/
void simulate(double dt);
void RunCells();
/**
* Returns the chemical field.
*/
auto GetField() const { return this->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)));
}
#ifndef POET_USE_PRM
struct SurrogateSetup {
std::vector<std::string> prop_names;
std::array<double, 2> base_totals;
bool has_het_ids;
/**
* 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);
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
@ -134,11 +121,42 @@ public:
* Enumerating DHT file options
*/
enum {
DHT_SNAPS_DISABLED = 0, //!< disabled file output
DHT_SNAPS_SIMEND, //!< only output of snapshot after simulation
DHT_SNAPS_ITEREND //!< output snapshots after each iteration
DHT_FILES_DISABLED = 0, //!< disabled file output
DHT_FILES_SIMEND, //!< only output of snapshot after simulation
DHT_FILES_ITEREND //!< output snapshots after each iteration
};
/**
* Initializes field with a "DataFrame" containing a single value for each
* species.
*
* Each species must have been previously found by
* ChemistryModule::RunInitFile.
*
* \exception std::domain_error Species name of map is not defined.
*
* \param n_cells Count of cells for each species.
* \param mapped_values Unordered map containing one double value for each
* specified species.
*/
void InitializeField(uint32_t n_cells, const SingleCMap &mapped_values);
/**
* Initializes field with a "DataFrame" containing a vector of values for each
* species.
*
* Each species must have been previously found by
* ChemistryModule::RunInitFile.
*
* There is no check if vector length matches count of grid cells defined.
*
* \exception std::domain_error Species name of map is not defined.
*
* \param mapped_values Unordered map containing a vector of multiple values.
* Size of the vectors shall be the count of grid cells defined previously.
*/
void InitializeField(const VectorCMap &mapped_values);
/**
* **Only called by workers!** Start the worker listening loop.
*/
@ -149,6 +167,55 @@ public:
*/
void MasterLoopBreak();
/**
* **Master only** Enables DHT usage for the workers.
*
* If called multiple times enabling the DHT, the current state of DHT will be
* overwritten with a new instance of the DHT.
*
* \param enable Enables or disables the usage of the DHT.
* \param size_mb Size in megabyte to allocate for the DHT if enabled.
*/
void SetDHTEnabled(bool enable, uint32_t size_mb);
/**
* **Master only** Set DHT snapshots to specific mode.
*
* \param type DHT snapshot mode.
* \param out_dir Path to output DHT snapshots.
*/
void SetDHTSnaps(int type, const std::string &out_dir);
/**
* **Master only** Set the vector with significant digits to round before
* inserting into DHT.
*
* \param signif_vec Vector defining significant digit for each species. Order
* is defined by prop_type vector (ChemistryModule::GetPropNames).
*/
void SetDHTSignifVector(std::vector<uint32_t> signif_vec);
/**
* **Master only** Set the DHT rounding type of each species. See
* DHT_PROP_TYPES enumemartion for explanation.
*
* \param proptype_vec Vector defining DHT prop type for each species.
*/
void SetDHTPropTypeVector(std::vector<uint32_t> proptype_vec);
/**
* **Master only** Load the state of the DHT from given file.
*
* \param input_file File to load data from.
*/
void ReadDHTFile(const std::string &input_file);
/**
* Overwrite the current field state by another field.
*
* There are no checks if new vector dimensions matches expected sizes etc.
*
* \param field New input field. Current field state of the instance will be
* overwritten.
*/
void SetField(const std::vector<double> &field) { this->field = field; }
/**
* **Master only** Return count of grid cells.
*/
@ -210,69 +277,35 @@ public:
*/
std::vector<uint32_t> GetWorkerDHTHits() const;
/**
* **Master only** Collect and return DHT misses of all workers.
*
* \return Vector of all count of DHT misses.
*/
std::vector<uint32_t> GetWorkerDHTMiss() const;
/**
* **Master only** Collect and return DHT evictions of all workers.
*
* \return Vector of all count of DHT evictions.
*/
std::vector<uint32_t> GetWorkerDHTEvictions() const;
/**
* **Master only** Returns the current state of the chemical field.
*
* \return Reference to the chemical field.
*/
Field &getField() { return this->chem_field; }
/**
* **Master only** Enable/disable progress bar.
*
* \param enabled True if print progressbar, false if not.
*/
void setProgressBarPrintout(bool enabled) {
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;
std::vector<double> GetWorkerInterpolationReadTimings() const;
std::vector<double> GetWorkerInterpolationGatherTimings() const;
std::vector<double> GetWorkerInterpolationFunctionCallTimings() const;
std::vector<uint32_t> GetWorkerPHTCacheHits() const;
std::vector<int> ai_surrogate_validity_vector;
#endif
protected:
void initializeDHT(uint32_t size_mb,
const NamedVector<std::uint32_t> &key_species,
bool has_het_ids);
void setDHTSnapshots(int type, const std::string &out_dir);
void setDHTReadFile(const std::string &input_file);
void initializeInterp(std::uint32_t bucket_size, std::uint32_t size_mb,
std::uint32_t min_entries,
const NamedVector<std::uint32_t> &key_species);
#ifdef POET_USE_PRM
void PrmToPoetField(std::vector<double> &field);
#else
enum {
CHEM_FIELD_INIT,
CHEM_INIT,
CHEM_WP_SIZE,
CHEM_DHT_ENABLE,
CHEM_DHT_SIGNIF_VEC,
CHEM_DHT_PROP_TYPE_VEC,
CHEM_DHT_SNAPS,
CHEM_DHT_READ_FILE,
CHEM_IP_ENABLE,
CHEM_IP_MIN_ENTRIES,
CHEM_IP_SIGNIF_VEC,
CHEM_WORK_LOOP,
CHEM_PERF,
CHEM_BREAK_MAIN_LOOP,
CHEM_AI_BCAST_VALIDITY
CHEM_BREAK_MAIN_LOOP
};
enum { LOOP_WORK, LOOP_END };
@ -282,20 +315,11 @@ protected:
WORKER_DHT_GET,
WORKER_DHT_FILL,
WORKER_IDLE,
WORKER_IP_WRITE,
WORKER_IP_READ,
WORKER_IP_GATHER,
WORKER_IP_FC,
WORKER_DHT_HITS,
WORKER_DHT_EVICTIONS,
WORKER_PHT_CACHE_HITS,
WORKER_IP_CALLS
WORKER_DHT_MISS,
WORKER_DHT_EVICTIONS
};
std::vector<uint32_t> interp_calls;
std::vector<uint32_t> dht_hits;
std::vector<uint32_t> dht_evictions;
struct worker_s {
double phreeqc_t = 0.;
double dht_get = 0.;
@ -311,7 +335,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,24 +361,16 @@ protected:
void WorkerPerfToMaster(int type, const struct worker_s &timings);
void WorkerMetricsToMaster(int type);
void WorkerRunWorkPackage(WorkPackage &work_package, double dSimTime,
double dTimestep);
IRM_RESULT WorkerRunWorkPackage(std::vector<double> &vecWP,
std::vector<int32_t> &vecMapping,
double dSimTime, double dTimestep);
void GetWPFromInternals(std::vector<double> &vecWP, uint32_t wp_size);
void SetInternalsFromWP(const std::vector<double> &vecWP, uint32_t wp_size);
std::vector<uint32_t> CalculateWPSizesVector(uint32_t n_cells,
uint32_t wp_size) const;
std::vector<double> shuffleField(const std::vector<double> &in_field,
uint32_t size_per_prop, uint32_t prop_count,
uint32_t wp_count);
void unshuffleField(const std::vector<double> &in_buffer,
uint32_t size_per_prop, uint32_t prop_count,
uint32_t wp_count, std::vector<double> &out_field);
std::vector<std::int32_t>
parseDHTSpeciesVec(const NamedVector<std::uint32_t> &key_species,
const std::vector<std::string> &to_compare) const;
void BCastStringVec(std::vector<std::string> &io);
int comm_size, comm_rank;
MPI_Comm group_comm;
@ -362,16 +378,11 @@ protected:
bool is_master;
uint32_t wp_size;
bool dht_enabled{false};
int dht_snaps_type{DHT_SNAPS_DISABLED};
bool dht_enabled = false;
int dht_snaps_type = DHT_FILES_DISABLED;
std::string dht_file_out_dir;
poet::DHT_Wrapper *dht = nullptr;
bool interp_enabled{false};
std::unique_ptr<poet::InterpolationModule> interp;
bool ai_surrogate_enabled{false};
poet::DHT_Wrapper *dht = std::nullptr_t();
static constexpr uint32_t BUFFER_OFFSET = 5;
@ -386,24 +397,16 @@ protected:
double idle_t = 0.;
double seq_t = 0.;
double send_recv_t = 0.;
#endif
std::array<double, 2> base_totals{0};
bool print_progessbar{false};
std::uint8_t file_pad{1};
double chem_t{0.};
double chem_t = 0.;
uint32_t n_cells = 0;
uint32_t prop_count = 0;
std::vector<std::string> prop_names;
Field chem_field;
const InitialList::ChemistryInit params;
std::unique_ptr<PhreeqcRunner> pqc_runner;
std::vector<double> field;
std::vector<uint32_t> speciesPerModule;
static constexpr uint32_t MODULE_COUNT = 5;
};
} // namespace poet

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

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

8
include/poet/DHT_Types.hpp Normal file
View File

@ -0,0 +1,8 @@
#ifndef DHT_TYPES_H_
#define DHT_TYPES_H_
namespace poet {
enum DHT_PROP_TYPES { DHT_TYPE_DEFAULT, DHT_TYPE_ACT, DHT_TYPE_IGNORE };
}
#endif // DHT_TYPES_H_

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

@ -1,4 +1,3 @@
/*
** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
** Potsdam)
@ -22,19 +21,8 @@
#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 "LookupKey.hpp"
#include <array>
#include <cstdint>
#include <limits>
#include <string>
#include <utility>
#include <vector>
extern "C" {
@ -45,7 +33,33 @@ extern "C" {
namespace poet {
using DHT_Location = std::pair<std::uint32_t, std::uint32_t>;
using DHT_Keyelement = struct keyelem {
std::int8_t exp : 8;
std::int64_t significant : 56;
};
using DHT_ResultObject = struct DHTResobj {
uint32_t length;
std::vector<std::vector<DHT_Keyelement>> keys;
std::vector<std::vector<double>> results;
std::vector<bool> needPhreeqc;
void ResultsToMapping(std::vector<int32_t> &curr_mapping);
void ResultsToWP(std::vector<double> &curr_wp);
};
/**
* @brief Return user-defined md5sum
*
* This function will calculate a hashsum with the help of md5sum. Therefore the
* md5sum for a given key is calculated and divided into two 64-bit parts. These
* will be XORed and returned as the hash.
*
* @param key_size Size of key in bytes
* @param key Pointer to key
* @return uint64_t Hashsum as an unsigned 64-bit integer
*/
static uint64_t get_md5(int key_size, void *key);
/**
* @brief C++-Wrapper around DHT implementation
@ -56,17 +70,6 @@ using DHT_Location = std::pair<std::uint32_t, std::uint32_t>;
*/
class DHT_Wrapper {
public:
using DHT_ResultObject = struct DHTResobj {
std::vector<LookupKey> keys;
std::vector<DHT_Location> locations;
std::vector<bool> filledDHT;
};
static constexpr std::int32_t DHT_KEY_INPUT_CUSTOM =
std::numeric_limits<std::int32_t>::min();
static constexpr int DHT_KEY_SIGNIF_DEFAULT = 5;
/**
* @brief Construct a new dht wrapper object
*
@ -76,18 +79,12 @@ public:
*
* @param dht_comm Communicator which addresses all participating DHT
* processes
* @param buckets_per_process Count of buckets to allocate for each
* process
* @param key_indices Vector indexing elements of one grid cell used
* for key creation.
* @param buckets_per_process Count of buckets to allocate for each process
* @param key_count Count of key entries
* @param data_count Count of data entries
*/
DHT_Wrapper(MPI_Comm dht_comm, std::uint64_t dht_size,
const NamedVector<std::uint32_t> &key_species,
const std::vector<std::int32_t> &key_indices,
const std::vector<std::string> &output_names,
const InitialList::ChemistryHookFunctions &hooks,
uint32_t data_count, bool with_interp, bool has_het_ids);
DHT_Wrapper(MPI_Comm dht_comm, uint32_t dht_size, uint32_t key_count,
uint32_t data_count);
/**
* @brief Destroy the dht wrapper object
*
@ -99,9 +96,6 @@ public:
*/
~DHT_Wrapper();
DHT_Wrapper &operator=(const DHT_Wrapper &) = delete;
DHT_Wrapper(const DHT_Wrapper &) = delete;
/**
* @brief Check if values of workpackage are stored in DHT
*
@ -118,7 +112,8 @@ public:
* @param[in,out] work_package Pointer to current work package
* @param dt Current timestep of simulation
*/
auto checkDHT(WorkPackage &work_package) -> const DHT_ResultObject &;
auto checkDHT(int length, double dt, const std::vector<double> &work_package)
-> poet::DHT_ResultObject;
/**
* @brief Write simulated values into DHT
@ -136,9 +131,8 @@ public:
* outputs of the PHREEQC simulation
* @param dt Current timestep of simulation
*/
void fillDHT(const WorkPackage &work_package);
void resultsToWP(std::vector<double> &work_package);
void fillDHT(int length, const DHT_ResultObject &dht_check_data,
const std::vector<double> &results);
/**
* @brief Dump current DHT state into file.
@ -178,6 +172,13 @@ public:
*/
auto getHits() { return this->dht_hits; };
/**
* @brief Get the Misses object
*
* @return uint64_t Count of read misses
*/
auto getMisses() { return this->dht_miss; };
/**
* @brief Get the Evictions object
*
@ -185,86 +186,27 @@ public:
*/
auto getEvictions() { return this->dht_evictions; };
void resetCounter() {
this->dht_hits = 0;
this->dht_evictions = 0;
}
void SetSignifVector(std::vector<uint32_t> signif_vec);
auto getDataCount() { return this->data_count; }
auto getCommunicator() { return this->communicator; }
DHT *getDHT() { return this->dht_object; };
DHT_ResultObject &getDHTResults() { return this->dht_results; }
const auto &getKeyElements() const { return this->input_key_elements; }
const auto &getKeySpecies() const { return this->key_species; }
void setBaseTotals(double tot_h, double tot_o) {
this->base_totals = {tot_h, tot_o};
}
std::uint32_t getInputCount() const {
return this->input_key_elements.size();
}
std::uint32_t getOutputCount() const { return this->data_count; }
inline void prepareKeys(const std::vector<std::vector<double>> &input_values,
double dt) {
dht_results.keys.resize(input_values.size());
for (std::size_t i = 0; i < input_values.size(); i++) {
if (this->hooks.dht_fuzz.isValid()) {
dht_results.keys[i] = fuzzForDHT_R(input_values[i], dt);
} else {
dht_results.keys[i] = fuzzForDHT(input_values[i], dt);
}
}
}
void SetPropTypeVector(std::vector<uint32_t> prop_type_vec);
private:
uint32_t key_count;
uint32_t data_count;
DHT *dht_object;
MPI_Comm communicator;
LookupKey fuzzForDHT(const std::vector<double> &cell, double dt);
LookupKey fuzzForDHT_R(const std::vector<double> &cell, double dt);
std::vector<double>
outputToInputAndRates(const std::vector<double> &old_results,
const std::vector<double> &new_results);
std::vector<double>
inputAndRatesToOutput(const std::vector<double> &dht_data,
const std::vector<double> &input_values);
std::vector<double> outputToRates(const std::vector<double> &old_results,
const std::vector<double> &new_results);
std::vector<double> ratesToOutput(const std::vector<double> &dht_data,
const std::vector<double> &input_values);
std::vector<DHT_Keyelement> fuzzForDHT(int var_count, void *key, double dt);
uint32_t dht_hits = 0;
uint32_t dht_miss = 0;
uint32_t dht_evictions = 0;
NamedVector<std::uint32_t> key_species;
std::vector<uint32_t> dht_signif_vector;
std::vector<uint32_t> dht_prop_type_vector;
std::vector<std::uint32_t> dht_signif_vector;
std::vector<std::uint32_t> dht_prop_type_vector;
std::vector<std::int32_t> input_key_elements;
const std::vector<std::string> &output_names;
const InitialList::ChemistryHookFunctions &hooks;
const bool with_interp;
DHT_ResultObject dht_results;
std::array<double, 2> base_totals{0};
bool has_het_ids{false};
static constexpr int DHT_KEY_SIGNIF_DEFAULT = 5;
static constexpr int DHT_KEY_SIGNIF_TOTALS = 12;
static constexpr int DHT_KEY_SIGNIF_CHARGE = 3;
};
} // namespace poet

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

@ -21,13 +21,17 @@
#ifndef DIFFUSION_MODULE_H
#define DIFFUSION_MODULE_H
#include "DataStructures/Field.hpp"
#include "Init/InitialList.hpp"
#include <sys/types.h>
#include "poet/SimParams.hpp"
#include <array>
#include <bits/stdint-uintn.h>
#include <cmath>
#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 +50,7 @@ public:
*
* @param R RRuntime object
*/
DiffusionModule(const InitialList::DiffusionInit &init_list, Field field)
: param_list(init_list), transport_field(field){};
DiffusionModule(poet::DiffusionParams diffu_args, Grid &grid_);
/**
* @brief Run simulation for one iteration
@ -57,6 +60,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
*
@ -64,28 +75,39 @@ public:
*/
double getTransportTime();
/**
* \brief Returns the current diffusion field.
*
* \return Reference to the diffusion field.
*/
Field &getField() { return this->transport_field; }
private:
/**
* @brief Instance of RRuntime
*
*/
// RRuntime &R;
InitialList::DiffusionInit param_list;
enum { DIM_1D = 1, DIM_2D };
Field transport_field;
void initialize(poet::DiffusionParams args);
void RoundToZero(double *field, uint32_t cell_count) const;
Grid &grid;
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;
poet::StateMemory *state;
uint32_t n_cells_per_prop;
/**
* @brief time spent for transport
*
*/
double transport_t = 0.;
double transport_t = 0.f;
};
} // namespace poet

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

@ -1,11 +1,11 @@
#ifndef DATASTRUCTURES_H_
#define DATASTRUCTURES_H_
#include <Rcpp.h>
#ifndef FIELD_H_
#define FIELD_H_
#include <cstddef>
#include <cstdint>
#include <iterator>
#include <string>
#include <unordered_map>
#include <vector>
namespace poet {
@ -23,21 +23,18 @@ using FieldColumn = std::vector<double>;
*/
class Field : private std::unordered_map<std::string, FieldColumn> {
public:
Field(){};
/**
* Creates a new instance of a field with fixed expected vector size.
*
* \param vec_s expected vector size of each component/column.
* \param vec_s expected vector size of each component/column. \
*/
Field(std::uint32_t vec_s) : req_vec_size(vec_s){};
Field(uint32_t vec_s) : req_vec_size(vec_s){};
/**
* Initializes instance with a 2D vector and according names for each columnn.
* There is no check if names were given multiple times. The order of name
* vector also defines the ordering of the output.
*
* \param vec_s expected vector size of each component/column.
* \param input 2D vector using STL semantic describing the current state of
* the field.
* \param prop_vec Name of each vector in the input. Shall match
@ -46,15 +43,14 @@ public:
* \exception std::runtime_error If prop_vec size doesn't match input vector
* size or column vectors size doesn't match expected vector size.
*/
Field(std::uint32_t vec_s, const std::vector<std::vector<double>> &input,
const std::vector<std::string> &prop_vec, bool row_major = false);
void InitFromVec(const std::vector<std::vector<double>> &input,
const std::vector<std::string> &prop_vec);
/**
* Initializes instance with a 1D continious memory vector and according names
* for each columnn. There is no check if names were given multiple times. The
* order of name vector also defines the ordering of the output.
*
* \param vec_s expected vector size of each component/column.
* \param input 1D vector using STL semantic describing the current state of
* the field storing each column starting at index *i times requested vector
* size*.
@ -64,52 +60,8 @@ public:
* \exception std::runtime_error If prop_vec size doesn't match input vector
* size or column vectors size doesn't match expected vector size.
*/
Field(std::uint32_t vec_s, const std::vector<double> &input,
const std::vector<std::string> &prop_vec);
Field(const SEXP &s_init) { fromSEXP(s_init); }
Field &operator=(const Field &rhs) {
this->req_vec_size = rhs.req_vec_size;
this->props = rhs.props;
this->clear();
for (const auto &pair : rhs) {
this->insert(pair);
}
return *this;
}
/**
* Read in a (previously exported) 1D vector. It has to have the same
* dimensions as the current column count times the requested vector size of
* this instance. Import order is given by the species name vector.
*
* \param cont_field 1D field as vector.
*
* \exception std::runtime_error Input vector does not match the expected
* size.
*/
Field &operator=(const FieldColumn &cont_field);
/**
* Read in a (previously exported) 2D vector. It has to have the same
* dimensions as the current column count of this instance and each vector
* must have the size of the requested vector size. Import order is given by
* the species name vector.
*
* \param cont_field 2D field as vector.
*
* \exception std::runtime_error Input vector has more or less elements than
* the instance or a column vector does not match expected vector size.
*/
Field &operator=(const std::vector<FieldColumn> &cont_field);
Field &operator=(const SEXP &s_rhs) {
fromSEXP(s_rhs);
return *this;
}
void InitFromVec(const std::vector<double> &input,
const std::vector<std::string> &prop_vec);
/**
* Returns a reference to the column vector with given name. Creates a new
@ -127,16 +79,14 @@ public:
*
* \return Vector containing all species names in output order.
*/
auto GetProps() const -> std::vector<std::string> { return this->props; }
auto GetProps() const { return this->props; };
/**
* Return the requested vector size.
*
* \return Requested vector size set in the instanciation of the object.
*/
auto GetRequestedVecSize() const -> std::uint32_t {
return this->req_vec_size;
};
auto GetRequestedVecSize() const { return this->req_vec_size; };
/**
* Updates all species with values from another field. If one element of the
@ -145,7 +95,7 @@ public:
*
* \param input Field to update the current instance's columns.
*/
void update(const Field &input);
void UpdateFromField(const Field &input);
/**
* Builds a new 1D vector with the current state of the instance. The output
@ -166,17 +116,36 @@ public:
*/
auto As2DVector() const -> std::vector<FieldColumn>;
SEXP asSEXP() const;
/**
* Read in a (previously exported) 1D vector. It has to have the same
* dimensions as the current column count times the requested vector size of
* this instance. Import order is given by the species name vector.
*
* \param cont_field 1D field as vector.
*
* \exception std::runtime_error Input vector does not match the expected
* size.
*/
void SetFromVector(const FieldColumn &cont_field);
std::size_t cols() const { return this->props.size(); }
/**
* Read in a (previously exported) 2D vector. It has to have the same
* dimensions as the current column count of this instance and each vector
* must have the size of the requested vector size. Import order is given by
* the species name vector.
*
* \param cont_field 2D field as vector.
*
* \exception std::runtime_error Input vector has more or less elements than
* the instance or a column vector does not match expected vector size.
*/
void SetFromVector(const std::vector<FieldColumn> &cont_field);
private:
std::uint32_t req_vec_size{0};
const uint32_t req_vec_size;
std::vector<std::string> props{{}};
void fromSEXP(const SEXP &s_rhs);
std::vector<std::string> props;
};
} // namespace poet
#endif // DATASTRUCTURES_H_
#endif // FIELD_H_

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

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

@ -1,4 +1,3 @@
/*
** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
** Potsdam)
@ -23,14 +22,17 @@
#define HASHFUNCTIONS_H_
#include <cstdint>
#include <openssl/evp.h>
namespace poet {
// Sum of POET interpreted as ASCII
constexpr uint32_t HASH_SEED = 80 + 79 + 69 + 84;
void initHashCtx(const EVP_MD *md);
void freeHashCtx();
uint64_t Murmur2_64A(int len, const void *key);
uint64_t hashDHT(int key_size, void *key);
} // namespace poet
}
#endif // HASHFUNCTIONS_H_

249
include/poet/SimParams.hpp Normal file
View File

@ -0,0 +1,249 @@
/*
** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
** Potsdam)
**
** Copyright (C) 2018-2022 Marco De Lucia, Max Luebke (GFZ Potsdam)
**
** POET is free software; you can redistribute it and/or modify it under the
** terms of the GNU General Public License as published by the Free Software
** Foundation; either version 2 of the License, or (at your option) any later
** version.
**
** POET is distributed in the hope that it will be useful, but WITHOUT ANY
** WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
** A PARTICULAR PURPOSE. See the GNU General Public License for more details.
**
** You should have received a copy of the GNU General Public License along with
** this program; if not, write to the Free Software Foundation, Inc., 51
** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
#ifndef PARSER_H
#define PARSER_H
#include <bits/stdint-uintn.h>
#include <string>
#include <string_view>
#include <vector>
#include "argh.hpp" // Argument handler https://github.com/adishavit/argh
#include <RInside.h>
#include <Rcpp.h>
// BSD-licenced
/** Standard DHT Size. Defaults to 1 GB (1000 MB) */
constexpr uint32_t DHT_SIZE_PER_PROCESS_MB = 1E3;
/** Standard work package size */
#define WORK_PACKAGE_SIZE_DEFAULT 5
namespace poet {
enum { PARSER_OK, PARSER_ERROR, PARSER_HELP };
/**
* @brief Defining all simulation parameters
*
*/
typedef struct {
/** Count of processes in MPI_COMM_WORLD */
int world_size;
/** rank of proces in MPI_COMM_WORLD */
int world_rank;
/** indicates if DHT should be used */
bool dht_enabled;
/** apply logarithm to key before rounding */
bool dht_log;
/** indicates if timestep dt differs between iterations */
bool dt_differ;
/** Indicates, when a DHT snapshot should be written */
int dht_snaps;
/** <b>not implemented</b>: How a DHT is distributed over processes */
int dht_strategy;
/** Size of DHt per process in byter */
unsigned int dht_size_per_process;
/** Default significant digit for rounding */
int dht_significant_digits;
/** Default work package size */
unsigned int wp_size;
/** indicates if resulting grid should be stored after every iteration */
bool store_result;
} t_simparams;
using GridParams = struct s_GridParams {
std::vector<uint32_t> n_cells;
std::vector<double> s_cells;
std::string type;
Rcpp::DataFrame init_df;
std::string input_script;
std::string database_path;
std::vector<std::string> props;
s_GridParams(RInside &R);
};
using DiffusionParams = struct s_DiffusionParams {
std::vector<std::string> prop_names;
Rcpp::NumericVector alpha;
Rcpp::List vecinj_inner;
Rcpp::DataFrame vecinj;
Rcpp::DataFrame vecinj_index;
s_DiffusionParams(RInside &R);
};
using ChemistryParams = struct s_ChemistryParams {
std::string database_path;
std::string input_script;
s_ChemistryParams(RInside &R);
};
/**
* @brief Reads information from program arguments and R runtime
*
* Providing functions to initialize parameters of the simulation using command
* line parameters and parameters from the R runtime. This class will also parse
* arguments from the commandline and decides if argument is known or unknown.
*
* Stores and distribute current simulation parameters at any time.
*
*/
class SimParams {
public:
/**
* @brief Construct a new SimParams object
*
* With all given parameters a new instance of this class will be created.
*
* @param world_rank Rank of process inside MPI_COMM_WORLD
* @param world_size Size of communicator MPI_COMM_WORLD
*/
SimParams(int world_rank, int world_size);
/**
* @brief Parse program arguments
*
* This is done by the argh.h library.
*
* First, the function will check if there is a flag 'help' or 'h'. If this is
* the case a help message is printed and the function will return with
* PARSER_HELP.
*
* Second, if there are not 2 positional arguments an error will be printed to
* stderr and the function returns with PARSER_ERROR.
*
* Then all given program parameters and flags will be read and checked, if
* there are known by validateOptions. A list of all unknown options might be
* returned, printed out and the function will return with PARSER_ERROR.
* Oterhwise the function continuos.
*
* Now all program arguments will be stored inside t_simparams struct, printed
* out and the function returns with PARSER_OK.
*
* Also, all parsed agruments are distributed to the R runtime.
*
* @param argv Argument value of the program
* @param R Instantiated R runtime
* @return int Returns with 0 if no error occured, otherwise value less than 0
* is returned.
*/
int parseFromCmdl(char *argv[], RInside &R);
/**
* @brief Init std::vector values
*
* This will initialize dht_signif_vector and dht_prop_type_vector internally
* depending on whether vectors are defined by R-Simulation file or not.
* 'init_chemistry' must be run beforehand.
*
* @param R R runtime
* @param col_count Count of variables per grid cell (typically the count of
* columns of each grid cell)
*/
void initVectorParams(RInside &R, int col_count);
/**
* @brief Get the numerical params struct
*
* Returns a struct which contains all numerical or boolean simulation
* parameters.
*
* @return t_simparams Parameter struct
*/
auto getNumParams() const { return this->simparams; };
/**
* @brief Get the DHT_Signif_Vector
*
* Returns a vector indicating which significant values are used for each
* variable of a grid cell.
*
* @return std::vector<int> Vector of integers containing information about
* significant digits
*/
auto getDHTSignifVector() const { return this->dht_signif_vector; };
/**
* @brief Get the DHT_Prop_Type_Vector
*
* Returns a vector indicating of which type a variable of a grid cell is.
*
* @return std::vector<std::string> Vector if strings defining a type of a
* variable
*/
auto getDHTPropTypeVector() const { return this->dht_prop_type_vector; };
/**
* @brief Return name of DHT snapshot.
*
* Name of the DHT file which is used to initialize the DHT with a previously
* written snapshot.
*
* @return std::string Absolute paht to the DHT snapshot
*/
auto getDHTFile() const { return this->dht_file; };
/**
* @brief Get the filesim name
*
* Returns a string containing the absolute path to a R file defining the
* simulation.
*
* @return std::string Absolute path to R file
*/
auto getFilesim() const { return this->filesim; };
/**
* @brief Get the output directory
*
* Returns the name of an absolute path where all output files should be
* stored.
*
* @return std::string Absolute path to output path
*/
auto getOutDir() const { return this->out_dir; };
private:
std::list<std::string> validateOptions(argh::parser cmdl);
const std::set<std::string> flaglist{"ignore-result", "dht", "dht-nolog"};
const std::set<std::string> paramlist{"work-package-size", "dht-signif",
"dht-strategy", "dht-size",
"dht-snaps", "dht-file"};
t_simparams simparams;
std::vector<uint32_t> dht_signif_vector;
std::vector<uint32_t> dht_prop_type_vector;
std::string dht_file;
std::string filesim;
std::string out_dir;
};
} // namespace poet
#endif // PARSER_H

459
include/poet/argh.hpp Normal file
View File

@ -0,0 +1,459 @@
/*
** Copyright (c) 2016, Adi Shavit All rights reserved.
**
** Redistribution and use in source and binary forms, with or without
** modification, are permitted provided that the following conditions are met:
**
** * Redistributions of source code must retain the above copyright notice, this
** list of conditions and the following disclaimer. * Redistributions in
** binary form must reproduce the above copyright notice, this list of
** conditions and the following disclaimer in the documentation and/or other
** materials provided with the distribution. * Neither the name of nor the
** names of its contributors may be used to endorse or promote products
** derived from this software without specific prior written permission.
**
** THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
** AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
** IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
** ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
** LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
** CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
** SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
** INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
** CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
** ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
** POSSIBILITY OF SUCH DAMAGE.
*/
#pragma once
#include <algorithm>
#include <sstream>
#include <string>
#include <vector>
#include <set>
#include <map>
#include <cassert>
namespace argh
{
// Terminology:
// A command line is composed of 2 types of args:
// 1. Positional args, i.e. free standing values
// 2. Options: args beginning with '-'. We identify two kinds:
// 2.1: Flags: boolean options => (exist ? true : false)
// 2.2: Parameters: a name followed by a non-option value
#if !defined(__GNUC__) || (__GNUC__ >= 5)
using string_stream = std::istringstream;
#else
// Until GCC 5, istringstream did not have a move constructor.
// stringstream_proxy is used instead, as a workaround.
class stringstream_proxy
{
public:
stringstream_proxy() = default;
// Construct with a value.
stringstream_proxy(std::string const& value) :
stream_(value)
{}
// Copy constructor.
stringstream_proxy(const stringstream_proxy& other) :
stream_(other.stream_.str())
{
stream_.setstate(other.stream_.rdstate());
}
void setstate(std::ios_base::iostate state) { stream_.setstate(state); }
// Stream out the value of the parameter.
// If the conversion was not possible, the stream will enter the fail state,
// and operator bool will return false.
template<typename T>
stringstream_proxy& operator >> (T& thing)
{
stream_ >> thing;
return *this;
}
// Get the string value.
std::string str() const { return stream_.str(); }
std::stringbuf* rdbuf() const { return stream_.rdbuf(); }
// Check the state of the stream.
// False when the most recent stream operation failed
operator bool() const { return !!stream_; }
~stringstream_proxy() = default;
private:
std::istringstream stream_;
};
using string_stream = stringstream_proxy;
#endif
class parser
{
public:
enum Mode { PREFER_FLAG_FOR_UNREG_OPTION = 1 << 0,
PREFER_PARAM_FOR_UNREG_OPTION = 1 << 1,
NO_SPLIT_ON_EQUALSIGN = 1 << 2,
SINGLE_DASH_IS_MULTIFLAG = 1 << 3,
};
parser() = default;
parser(std::initializer_list<char const* const> pre_reg_names)
{ add_params(pre_reg_names); }
parser(const char* const argv[], int mode = PREFER_FLAG_FOR_UNREG_OPTION)
{ parse(argv, mode); }
parser(int argc, const char* const argv[], int mode = PREFER_FLAG_FOR_UNREG_OPTION)
{ parse(argc, argv, mode); }
void add_param(std::string const& name);
void add_params(std::initializer_list<char const* const> init_list);
void parse(const char* const argv[], int mode = PREFER_FLAG_FOR_UNREG_OPTION);
void parse(int argc, const char* const argv[], int mode = PREFER_FLAG_FOR_UNREG_OPTION);
std::multiset<std::string> const& flags() const { return flags_; }
std::map<std::string, std::string> const& params() const { return params_; }
std::vector<std::string> const& pos_args() const { return pos_args_; }
// begin() and end() for using range-for over positional args.
std::vector<std::string>::const_iterator begin() const { return pos_args_.cbegin(); }
std::vector<std::string>::const_iterator end() const { return pos_args_.cend(); }
size_t size() const { return pos_args_.size(); }
//////////////////////////////////////////////////////////////////////////
// Accessors
// flag (boolean) accessors: return true if the flag appeared, otherwise false.
bool operator[](std::string const& name) const;
// multiple flag (boolean) accessors: return true if at least one of the flag appeared, otherwise false.
bool operator[](std::initializer_list<char const* const> init_list) const;
// returns positional arg string by order. Like argv[] but without the options
std::string const& operator[](size_t ind) const;
// returns a std::istream that can be used to convert a positional arg to a typed value.
string_stream operator()(size_t ind) const;
// same as above, but with a default value in case the arg is missing (index out of range).
template<typename T>
string_stream operator()(size_t ind, T&& def_val) const;
// parameter accessors, give a name get an std::istream that can be used to convert to a typed value.
// call .str() on result to get as string
string_stream operator()(std::string const& name) const;
// accessor for a parameter with multiple names, give a list of names, get an std::istream that can be used to convert to a typed value.
// call .str() on result to get as string
// returns the first value in the list to be found.
string_stream operator()(std::initializer_list<char const* const> init_list) const;
// same as above, but with a default value in case the param was missing.
// Non-string def_val types must have an operator<<() (output stream operator)
// If T only has an input stream operator, pass the string version of the type as in "3" instead of 3.
template<typename T>
string_stream operator()(std::string const& name, T&& def_val) const;
// same as above but for a list of names. returns the first value to be found.
template<typename T>
string_stream operator()(std::initializer_list<char const* const> init_list, T&& def_val) const;
private:
string_stream bad_stream() const;
std::string trim_leading_dashes(std::string const& name) const;
bool is_number(std::string const& arg) const;
bool is_option(std::string const& arg) const;
bool got_flag(std::string const& name) const;
bool is_param(std::string const& name) const;
private:
std::vector<std::string> args_;
std::map<std::string, std::string> params_;
std::vector<std::string> pos_args_;
std::multiset<std::string> flags_;
std::set<std::string> registeredParams_;
std::string empty_;
};
//////////////////////////////////////////////////////////////////////////
inline void parser::parse(const char * const argv[], int mode)
{
int argc = 0;
for (auto argvp = argv; *argvp; ++argc, ++argvp);
parse(argc, argv, mode);
}
//////////////////////////////////////////////////////////////////////////
inline void parser::parse(int argc, const char* const argv[], int mode /*= PREFER_FLAG_FOR_UNREG_OPTION*/)
{
// convert to strings
args_.resize(argc);
std::transform(argv, argv + argc, args_.begin(), [](const char* const arg) { return arg; });
// parse line
for (auto i = 0u; i < args_.size(); ++i)
{
if (!is_option(args_[i]))
{
pos_args_.emplace_back(args_[i]);
continue;
}
auto name = trim_leading_dashes(args_[i]);
if (!(mode & NO_SPLIT_ON_EQUALSIGN))
{
auto equalPos = name.find('=');
if (equalPos != std::string::npos)
{
params_.insert({ name.substr(0, equalPos), name.substr(equalPos + 1) });
continue;
}
}
// if the option is unregistered and should be a multi-flag
if (1 == (args_[i].size() - name.size()) && // single dash
argh::parser::SINGLE_DASH_IS_MULTIFLAG & mode && // multi-flag mode
!is_param(name)) // unregistered
{
std::string keep_param;
if (!name.empty() && is_param(std::string(1ul, name.back()))) // last char is param
{
keep_param += name.back();
name.resize(name.size() - 1);
}
for (auto const& c : name)
{
flags_.emplace(std::string{ c });
}
if (!keep_param.empty())
{
name = keep_param;
}
else
{
continue; // do not consider other options for this arg
}
}
// any potential option will get as its value the next arg, unless that arg is an option too
// in that case it will be determined a flag.
if (i == args_.size() - 1 || is_option(args_[i + 1]))
{
flags_.emplace(name);
continue;
}
// if 'name' is a pre-registered option, then the next arg cannot be a free parameter to it is skipped
// otherwise we have 2 modes:
// PREFER_FLAG_FOR_UNREG_OPTION: a non-registered 'name' is determined a flag.
// The following value (the next arg) will be a free parameter.
//
// PREFER_PARAM_FOR_UNREG_OPTION: a non-registered 'name' is determined a parameter, the next arg
// will be the value of that option.
assert(!(mode & argh::parser::PREFER_FLAG_FOR_UNREG_OPTION)
|| !(mode & argh::parser::PREFER_PARAM_FOR_UNREG_OPTION));
bool preferParam = mode & argh::parser::PREFER_PARAM_FOR_UNREG_OPTION;
if (is_param(name) || preferParam)
{
params_.insert({ name, args_[i + 1] });
++i; // skip next value, it is not a free parameter
continue;
}
else
{
flags_.emplace(name);
}
};
}
//////////////////////////////////////////////////////////////////////////
inline string_stream parser::bad_stream() const
{
string_stream bad;
bad.setstate(std::ios_base::failbit);
return bad;
}
//////////////////////////////////////////////////////////////////////////
inline bool parser::is_number(std::string const& arg) const
{
// inefficient but simple way to determine if a string is a number (which can start with a '-')
std::istringstream istr(arg);
double number;
istr >> number;
return !(istr.fail() || istr.bad());
}
//////////////////////////////////////////////////////////////////////////
inline bool parser::is_option(std::string const& arg) const
{
assert(0 != arg.size());
if (is_number(arg))
return false;
return '-' == arg[0];
}
//////////////////////////////////////////////////////////////////////////
inline std::string parser::trim_leading_dashes(std::string const& name) const
{
auto pos = name.find_first_not_of('-');
return std::string::npos != pos ? name.substr(pos) : name;
}
//////////////////////////////////////////////////////////////////////////
inline bool argh::parser::got_flag(std::string const& name) const
{
return flags_.end() != flags_.find(trim_leading_dashes(name));
}
//////////////////////////////////////////////////////////////////////////
inline bool argh::parser::is_param(std::string const& name) const
{
return registeredParams_.count(name);
}
//////////////////////////////////////////////////////////////////////////
inline bool parser::operator[](std::string const& name) const
{
return got_flag(name);
}
//////////////////////////////////////////////////////////////////////////
inline bool parser::operator[](std::initializer_list<char const* const> init_list) const
{
return std::any_of(init_list.begin(), init_list.end(), [&](char const* const name) { return got_flag(name); });
}
//////////////////////////////////////////////////////////////////////////
inline std::string const& parser::operator[](size_t ind) const
{
if (ind < pos_args_.size())
return pos_args_[ind];
return empty_;
}
//////////////////////////////////////////////////////////////////////////
inline string_stream parser::operator()(std::string const& name) const
{
auto optIt = params_.find(trim_leading_dashes(name));
if (params_.end() != optIt)
return string_stream(optIt->second);
return bad_stream();
}
//////////////////////////////////////////////////////////////////////////
inline string_stream parser::operator()(std::initializer_list<char const* const> init_list) const
{
for (auto& name : init_list)
{
auto optIt = params_.find(trim_leading_dashes(name));
if (params_.end() != optIt)
return string_stream(optIt->second);
}
return bad_stream();
}
//////////////////////////////////////////////////////////////////////////
template<typename T>
string_stream parser::operator()(std::string const& name, T&& def_val) const
{
auto optIt = params_.find(trim_leading_dashes(name));
if (params_.end() != optIt)
return string_stream(optIt->second);
std::ostringstream ostr;
ostr << def_val;
return string_stream(ostr.str()); // use default
}
//////////////////////////////////////////////////////////////////////////
// same as above but for a list of names. returns the first value to be found.
template<typename T>
string_stream parser::operator()(std::initializer_list<char const* const> init_list, T&& def_val) const
{
for (auto& name : init_list)
{
auto optIt = params_.find(trim_leading_dashes(name));
if (params_.end() != optIt)
return string_stream(optIt->second);
}
std::ostringstream ostr;
ostr << def_val;
return string_stream(ostr.str()); // use default
}
//////////////////////////////////////////////////////////////////////////
inline string_stream parser::operator()(size_t ind) const
{
if (pos_args_.size() <= ind)
return bad_stream();
return string_stream(pos_args_[ind]);
}
//////////////////////////////////////////////////////////////////////////
template<typename T>
string_stream parser::operator()(size_t ind, T&& def_val) const
{
if (pos_args_.size() <= ind)
{
std::ostringstream ostr;
ostr << def_val;
return string_stream(ostr.str());
}
return string_stream(pos_args_[ind]);
}
//////////////////////////////////////////////////////////////////////////
inline void parser::add_param(std::string const& name)
{
registeredParams_.insert(trim_leading_dashes(name));
}
//////////////////////////////////////////////////////////////////////////
inline void parser::add_params(std::initializer_list<char const* const> init_list)
{
for (auto& name : init_list)
registeredParams_.insert(trim_leading_dashes(name));
}
}

16
src/Base/Macros.hpp
View File

@ -1,16 +0,0 @@
#ifndef MACROS_H
#define MACROS_H
#include <iostream>
#include <string>
// Prepend "msg" with name of calling function
#define MSG(msg) std::cout << "CPP: " << __func__ << ": " << std::string(msg) << std::endl;
#define MSG_NOENDL(msg) std::cout << "CPP: " << __func__ << ": " << std::string(msg);
#define ERRMSG(msg) std::cerr << "CPP_ERROR: " << __func__ << ": " << std::string(msg) << std::endl;
#define BOOL_PRINT(bool) std::string(bool ? "ON" : "OFF")
#endif // MACROS_H

92
src/Base/RInsidePOET.hpp
View File

@ -1,92 +0,0 @@
#pragma once
#include <RInside.h>
#include <Rcpp.h>
#include <Rinternals.h>
#include <exception>
#include <memory>
#include <string>
namespace poet {
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(){};
};
/**
* @brief Deferred evaluation function
*
* The class is intended to call R functions within an existing RInside
* instance. The problem with "original" Rcpp::Function is that they require:
* 1. RInside instance already present, restricting the declaration of
* Rcpp::Functions in global scope
* 2. Require the function to be present. Otherwise, they will throw an
* exception.
* This class solves both problems by deferring the evaluation of the function
* until the constructor is called and evaluating whether the function is
* present or not, wihout throwing an exception.
*
* @tparam T Return type of the function
*/
class DEFunc {
public:
DEFunc() {}
DEFunc(const std::string &f_name) {
try {
this->func = std::make_shared<Rcpp::Function>(f_name);
} catch (const std::exception &e) {
}
}
DEFunc(SEXP f) {
try {
this->func = std::make_shared<Rcpp::Function>(f);
} catch (const std::exception &e) {
}
}
template <typename... Args> SEXP operator()(Args... args) const {
if (func) {
return (*this->func)(args...);
} else {
throw std::exception();
}
}
DEFunc &operator=(const DEFunc &rhs) {
this->func = rhs.func;
return *this;
}
DEFunc(const DEFunc &rhs) { this->func = rhs.func; }
bool isValid() const { return static_cast<bool>(func); }
SEXP asSEXP() const {
if (!func) {
return R_NilValue;
}
return Rcpp::as<SEXP>(*this->func.get());
}
private:
std::shared_ptr<Rcpp::Function> func;
};
} // namespace poet

109
src/CMakeLists.txt
View File

@ -1,101 +1,16 @@
add_library(POETLib
Init/InitialList.cpp
Init/GridInit.cpp
Init/DiffusionInit.cpp
Init/ChemistryInit.cpp
DataStructures/Field.cpp
Transport/DiffusionModule.cpp
Chemistry/ChemistryModule.cpp
Chemistry/MasterFunctions.cpp
Chemistry/WorkerFunctions.cpp
Chemistry/SurrogateModels/DHT_Wrapper.cpp
Chemistry/SurrogateModels/DHT.c
Chemistry/SurrogateModels/HashFunctions.cpp
Chemistry/SurrogateModels/InterpolationModule.cpp
Chemistry/SurrogateModels/ProximityHashTable.cpp
)
add_subdirectory(DataStructures)
set(POET_TUG_APPROACH "Implicit" CACHE STRING "tug numerical approach to use")
set_property(CACHE POET_TUG_APPROACH PROPERTY STRINGS "Implicit" "Explicit")
file(GLOB poet_lib_SRC
CONFIGURE_DEPENDS
"*.cpp" "*.c")
if (POET_TUG_APPROACH STREQUAL "Implicit")
target_compile_definitions(POETLib PRIVATE POET_TUG_BTCS)
elseif (POET_TUG_APPROACH STREQUAL "Explicit")
target_compile_definitions(POETLib PRIVATE POET_TUG_FTCS)
endif()
find_library(MATH_LIBRARY m)
find_library(CRYPTO_LIBRARY crypto)
target_include_directories(POETLib PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}")
target_link_libraries(
POETLib
PUBLIC RRuntime
PUBLIC IPhreeqcPOET
PUBLIC tug
PUBLIC MPI::MPI_C
)
add_library(poet_lib ${poet_lib_SRC})
target_include_directories(poet_lib PUBLIC ${PROJECT_SOURCE_DIR}/include)
target_link_libraries(poet_lib PUBLIC
MPI::MPI_CXX ${MATH_LIBRARY} ${CRYPTO_LIBRARY} RRuntime tug PhreeqcRM DataStructures ChemistryModule)
target_compile_definitions(poet_lib PUBLIC STRICT_R_HEADERS OMPI_SKIP_MPICXX)
include(FetchContent)
FetchContent_Declare(
cli11
QUIET
GIT_REPOSITORY https://github.com/CLIUtils/CLI11.git
GIT_TAG v2.5.0
)
FetchContent_MakeAvailable(cli11)
# add_library(poetlib
# Base/Grid.cpp
# Base/SimParams.cpp
# Chemistry/ChemistryModule.cpp
# Chemistry/MasterFunctions.cpp
# Chemistry/WorkerFunctions.cpp
# Chemistry/SurrogateModels/DHT_Wrapper.cpp
# Chemistry/SurrogateModels/DHT.c
# Chemistry/SurrogateModels/HashFunctions.cpp
# Chemistry/SurrogateModels/InterpolationModule.cpp
# Chemistry/SurrogateModels/ProximityHashTable.cpp
# Transport/DiffusionModule.cpp
# )
# target_link_libraries(poetlib PUBLIC
# DataStructures
# Init
# MPI::MPI_C
# ${MATH_LIBRARY}
# RRuntime
# PhreeqcRM
# tug
# )
target_compile_definitions(POETLib PUBLIC STRICT_R_HEADERS OMPI_SKIP_MPICXX)
# mark_as_advanced(PHREEQCRM_BUILD_MPI PHREEQCRM_DISABLE_OPENMP)
# set(PHREEQCRM_DISABLE_OPENMP ON CACHE BOOL "" FORCE)
# option(POET_DHT_DEBUG "Build with DHT debug info" OFF)
# if(POET_DHT_DEBUG)
# target_compile_definitions(poetlib PRIVATE DHT_STATISTICS)
# endif()
# option(POET_PHT_ADDITIONAL_INFO "Enables additional information in the PHT" OFF)
# if (POET_PHT_ADDITIONAL_INFO)
# target_compile_definitions(poetlib PRIVATE POET_PHT_ADD)
# endif()
file(READ "${PROJECT_SOURCE_DIR}/R_lib/kin_r_library.R" R_KIN_LIB )
file(READ "${PROJECT_SOURCE_DIR}/R_lib/init_r_lib.R" R_INIT_LIB)
file(READ "${PROJECT_SOURCE_DIR}/R_lib/ai_surrogate_model.R" R_AI_SURROGATE_LIB)
configure_file(poet.hpp.in poet.hpp @ONLY)
add_executable(poet poet.cpp)
target_link_libraries(poet PRIVATE POETLib MPI::MPI_C RRuntime CLI11::CLI11)
target_include_directories(poet PRIVATE "${CMAKE_CURRENT_BINARY_DIR}")
add_executable(poet_init initializer.cpp)
target_link_libraries(poet_init PRIVATE POETLib RRuntime CLI11::CLI11)
target_include_directories(poet_init PRIVATE "${CMAKE_CURRENT_BINARY_DIR}")
install(TARGETS poet poet_init DESTINATION bin)
add_subdirectory(ChemistryModule)

20
src/Chemistry/ChemistryDefs.hpp
View File

@ -1,20 +0,0 @@
#pragma once
#include <cstdint>
#include <vector>
namespace poet {
enum DHT_PROP_TYPES { DHT_TYPE_DEFAULT, DHT_TYPE_CHARGE, DHT_TYPE_TOTAL };
enum CHEMISTRY_OUT_SOURCE { CHEM_PQC, CHEM_DHT, CHEM_INTERP, CHEM_AISURR };
struct WorkPackage {
std::size_t size;
std::vector<std::vector<double>> input;
std::vector<std::vector<double>> output;
std::vector<std::uint8_t> mapping;
WorkPackage(std::size_t _size)
: size(_size), input(size), output(size), mapping(size, CHEM_PQC) {}
};
} // namespace poet

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