Merge branch 'mdl/defaultstorelease' into 'main'

DEFAULT_BUILD_TYPE set to Release

See merge request naaice/poet!60

.devcontainer/Dockerfile

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

.devcontainer/devcontainer.json

@@ -0,0 +1,29 @@
+// 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"
+}

.gitignore

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

.gitlab-ci.yml

@@ -19,7 +19,6 @@
 image: git.gfz-potsdam.de:5000/naaice/poet:ci
 
 stages:          # List of stages for jobs, and their order of execution
-  - build
   - release
   - test
 
@@ -28,89 +27,93 @@ variables:
   SOURCE_ARCHIVE_NAME: 'poet_${CI_COMMIT_TAG}_sources.tar.gz'
   CHANGELOG_FILE: 'commit_changelog.md'
 
-build-poet:       # This job runs in the build stage, which runs first.
-  stage: build
-  script:
-    - mkdir build && cd build
-    - cmake -DPOET_ENABLE_TESTING=ON ..
-    - make -j$(nproc)
-  artifacts:
-    paths:
-      - build
-    expire_in: 1 day
-
-test-poet:
+test:       # This job runs in the build stage, which runs first.
   stage: test
-  dependencies:
-    - build-poet
   script:
-    - cd build
+    - mkdir -p build && cd build
+    - cmake -DPOET_ENABLE_TESTING=ON -DPOET_PREPROCESS_BENCHS=OFF -DCMAKE_BUILD_TYPE=Release ..
     - make -j$(nproc) check
 
-archive-sources:       # This job runs in the build stage, which runs first.
-  image: python:3
-  stage: release
-
-  before_script:
-    - pip install git-archive-all
-    - echo ARCHIVE_JOB_ID=${CI_JOB_ID} >> archives.env
-  script:
-    - git-archive-all ${SOURCE_ARCHIVE_NAME}
-  artifacts:
-    paths:
-      - ${SOURCE_ARCHIVE_NAME}
-    expire_in: never
-    reports:
-      dotenv: archives.env
-  rules:
-    - if: $CI_COMMIT_TAG
-
-release-description:
-  image: golang:bullseye
-  stage: release
-  rules:
-    - if: $CI_COMMIT_TAG
-  before_script:
-    - go install github.com/git-chglog/git-chglog/cmd/git-chglog@v0.15.2
-  script:
-    - git-chglog -o ${CHANGELOG_FILE} ${CI_COMMIT_TAG}
-  artifacts:
-    paths:
-      - ${CHANGELOG_FILE}
-
-
-release-create:
-  stage: release
-  image: registry.gitlab.com/gitlab-org/release-cli:latest
-  rules:
-    - if: $CI_COMMIT_TAG
-  script:
-    - echo "Running release job"
-  needs:
-    - job: archive-sources
-      artifacts: true
-    - job: release-description
-      artifacts: true
-  release:
-    tag_name: $CI_COMMIT_TAG
-    name: 'POET $CI_COMMIT_TAG'
-    description: ${CHANGELOG_FILE}
-    assets:
-      links:
-        - name: '${SOURCE_ARCHIVE_NAME}'
-          url: 'https://git.gfz-potsdam.de/naaice/poet/-/jobs/${ARCHIVE_JOB_ID}/artifacts/file/${SOURCE_ARCHIVE_NAME}'
-
 pages:
   stage: release
   before_script:
     - apt-get update && apt-get install -y doxygen graphviz
-    - mkdir {build,public}
+    - mkdir {build_pages,public}
   script:
-    - pushd build
+    - pushd build_pages
     - cmake .. && make doxygen
-    - popd && mv build/docs/html/* public/
+    - popd && mv build_pages/docs/html/* public/
   artifacts:
     paths:
       - public
   rules:
-    - if: $CI_COMMIT_REF_NAME == $CI_DEFAULT_BRANCH
+    - 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}'

.gitmodules

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

CITATION.cff

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

CMake/FindRRuntime.cmake

@@ -1,30 +1,27 @@
 # prepare R environment (Rcpp + RInside)
-find_program(R_EXE "R")
+find_program(R_EXE "R"
+  REQUIRED
+)
 
 # search for R executable, R header file and library path
-if(R_EXE)
-  execute_process(
-    COMMAND ${R_EXE} RHOME
-    OUTPUT_VARIABLE R_ROOT_DIR
-    OUTPUT_STRIP_TRAILING_WHITESPACE
-  )
+execute_process(
+  COMMAND ${R_EXE} RHOME
+  OUTPUT_VARIABLE R_ROOT_DIR
+  OUTPUT_STRIP_TRAILING_WHITESPACE
+)
 
-  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_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_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)
+find_library(
+  R_LIBRARY libR.so
+  HINTS ${R_ROOT_DIR}/lib
+  REQUIRED
+)
 
 set(R_LIBRARIES ${R_LIBRARY})
 set(R_INCLUDE_DIRS ${R_INCLUDE_DIR})
@@ -45,8 +42,6 @@ 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
@@ -72,16 +67,11 @@ 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)
-set_target_properties(
-  RRuntime PROPERTIES
-  INTERFACE_LINK_LIBRARIES "${R_LIBRARIES}"
-  INTERFACE_INCLUDE_DIRECTORIES "${R_INCLUDE_DIRS}"
-)
+target_link_libraries(RRuntime INTERFACE ${R_LIBRARIES})
+target_include_directories(RRuntime INTERFACE ${R_INCLUDE_DIRS})
 
 unset(R_LIBRARIES)
 unset(R_INCLUDE_DIRS)

CMake/POET_Scripts.cmake

@@ -13,7 +13,7 @@ macro(get_POET_version)
       WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}
       OUTPUT_VARIABLE POET_GIT_VERSION
       OUTPUT_STRIP_TRAILING_WHITESPACE)
-    if(NOT POET_GIT_BRANCH STREQUAL "master")
+    if(NOT POET_GIT_BRANCH STREQUAL "main")
       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 "0.1")
+    set(POET_VERSION "not_versioned")
   endif()
 
   message(STATUS "Configuring POET version ${POET_VERSION}")

CMakeLists.txt

@@ -1,14 +1,23 @@
-# Version 3.9+ offers new MPI package variables
-cmake_minimum_required(VERSION 3.9)
+cmake_minimum_required(VERSION 3.20)
 
 project(POET
   LANGUAGES CXX C
   DESCRIPTION "A coupled reactive transport simulator")
 
 # specify the C++ standard
-set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD 20)
 set(CMAKE_CXX_STANDARD_REQUIRED True)
 
+set(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")
@@ -16,27 +25,27 @@ list(APPEND CMAKE_MODULE_PATH "${POET_SOURCE_DIR}/CMake")
 
 get_poet_version()
 
-# set(GCC_CXX_FLAGS "-D STRICT_R_HEADERS") add_definitions(${GCC_CXX_FLAGS})
-
 find_package(MPI REQUIRED)
 
 find_package(RRuntime REQUIRED)
 
 add_subdirectory(src)
-add_subdirectory(R_lib)
-add_subdirectory(app)
-add_subdirectory(bench)
+
+option(POET_PREPROCESS_BENCHS "Preprocess benchmarks" ON)
+
+if (POET_PREPROCESS_BENCHS)
+  add_subdirectory(bench)
+endif()
 
 # 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/phreeqcrm EXCLUDE_FROM_ALL)
+add_subdirectory(ext/litephreeqc 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()
 

LICENSE

@@ -1,8 +1,8 @@
                     GNU GENERAL PUBLIC LICENSE
                        Version 2, June 1991
 
- Copyright (C) 1989, 1991 Free Software Foundation, Inc., <http://fsf.org/>
- 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ Copyright (C) 1989, 1991 Free Software Foundation, Inc.,
+ <https://fsf.org/>
  Everyone is permitted to copy and distribute verbatim copies
  of this license document, but changing it is not allowed.
 
@@ -278,3 +278,61 @@ 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.

README.md

@@ -1,77 +1,108 @@
-<!--
-    Time-stamp: "Last modified 2023-08-02 13:55:11 mluebke"
--->
-
 # POET
 
-[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.
+**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's Coupling Scheme](./docs/20230720_Scheme_POET_en.svg)
+[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)
 
 ## Parsed code documentiation
 
-A parsed version of POET's documentiation can be found at [Gitlab
-pages](https://naaice.git-pages.gfz-potsdam.de/poet).
+A parsed version of POET's documentation can be found at [Gitlab
+pages](https://naaice.git-pages.gfz.de/poet).
 
 ## External Libraries
 
-The following external header library is shipped with POET:
+The following external libraries are shipped with POET:
 
-- **argh** - https://github.com/adishavit/argh (BSD license)
-- **PhreeqcRM** with patches from GFZ -
-  https://www.usgs.gov/software/phreeqc-version-3 -
-  https://git.gfz-potsdam.de/mluebke/phreeqcrm-gfz
-- **tug** - https://git.gfz-potsdam.de/sec34/tug
+- **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>
 
 ## Installation
 
 ### Requirements
 
-To compile POET you need several software to be installed:
+To compile POET you need following software to be installed:
 
 - C/C++ compiler (tested with GCC)
 - MPI-Implementation (tested with OpenMPI and MVAPICH)
-- R language and environment
-- CMake 3.9+
-- *optional*: `doxygen` with `dot` bindings for documentiation
+- 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)
 
-The following R libraries must then be installed, which will get the
-needed dependencies automatically:
+The following R packages (and their dependencies) must also be
+installed:
 
 - [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
 
-The generation of makefiles is done with CMake. You should be able to generate
-Makefiles by running:
+POET is built with CMake. You can generate Makefiles by running the
+usual:
 
 ```sh
 mkdir build && cd build
-cmake ..
+cmake -DCMAKE_BUILD_TYPE=Release ..
 ```
 
 This will create the directory `build` and processes the CMake files
 and generate Makefiles from it. You're now able to run `make` to start
 build process.
 
-If everything went well you'll find the executable at
-`build/app/poet`, but it is recommended to install the POET project
+If everything went well you'll find the executables at
+`build/src/poet`, but it is recommended to install the POET project
 structure to a desired `CMAKE_INSTALL_PREFIX` with `make install`.
 
 During the generation of Makefiles, various options can be specified
 via `cmake -D <option>=<value> [...]`. Currently, there are the
 following available options:
 
-- **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_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*.
   
 ### Example: Build from scratch
 
@@ -84,7 +115,7 @@ follows:
 $ R
 
 # install R dependencies
-> install.packages(c("Rcpp", "RInside"))
+> install.packages(c("Rcpp", "RInside","qs","qs2"))
 > q(save="no")
 
 # cd into POET project root
@@ -92,7 +123,7 @@ $ cd <POET_dir>
 
 # Build process
 $ mkdir build && cd build
-$ cmake -DCMAKE_INSTALL_PREFIX=/home/<user>/poet ..
+$ cmake -DCMAKE_INSTALL_PREFIX=/home/<user>/poet -DCMAKE_BUILD_TYPE=Release ..
 $ make -j<max_numprocs>
 $ make install
 ```
@@ -107,66 +138,64 @@ The correspondending directory tree would look like this:
 ```sh
 poet
 ├── bin
-│   └── poet
-├── R_lib
-│   └── kin_r_library.R
+│   ├── poet
+│   └── poet_init
 └── share
     └── poet
-        └── bench
-            ├── barite
-            │   ├── barite_interp_eval.R
-            │   ├── barite.pqi
-            │   ├── barite.R
-            │   └── db_barite.dat
-            ├── dolo
-            │   ├── dolo_diffu_inner_large.R
-            │   ├── dolo_diffu_inner.R
-            │   ├── dolo_inner.pqi
-            │   ├── dolo_interp_long.R
-            │   └── phreeqc_kin.dat
-            └── surfex
-                ├── ExBase.pqi
-                ├── ex.R
-                ├── SMILE_2021_11_01_TH.dat
-                ├── SurfExBase.pqi
-                └── surfex.R
+        ├── barite
+        │   ├── barite_200.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
 ```
 
-The R libraries will be loaded at runtime and the paths are hardcoded
-absolute paths inside `poet.cpp`. So, if you consider to move
-`bin/poet` either change paths of the R source files and recompile
-POET or also move `R_lib/*` relative to the binary.
+With the installation of POET, two executables are provided: 
+  - `poet` - the main executable to run simulations
+  - `poet_init` - a preprocessor to generate input files for POET from
+    R scripts
 
-The benchmarks consist of input scripts, which are provided as .R files.
-Additionally, Phreeqc scripts and their corresponding databases are required,
-stored as .pqi and .dat files, respectively.
+Preprocessed benchmarks can be found in the `share/poet` directory
+with an according *runtime* setup. More on those files and how to
+create them later.
 
 ## Running
 
-Run POET by `mpirun ./poet <OPTIONS> <SIMFILE> <OUTPUT_DIRECTORY>`
-where:
+Run POET by `mpirun ./poet [OPTIONS] <RUNFILE> <SIMFILE>
+<OUTPUT_DIRECTORY>` where:
 
-- **OPTIONS** - runtime parameters (explained below)
-- **SIMFILE** - simulation described as R script (e.g.
-  `<POET_INSTALL_DIR>/share/poet/bench/dolo/dolo_interp_long.R`)
-- **OUTPUT_DIRECTORY** - path, where all output of POET should be stored
+- **OPTIONS** - POET options (explained below)
+- **RUNFILE** - Runtime parameters described as R script
+- **SIMFILE** - Simulation input prepared by `poet_init`
+- **OUTPUT_DIRECTORY** - path, where all output of POET should be
+  stored
 
-### Runtime options
+### POET command line arguments
 
 The following parameters can be set:
 
-| Option                      | Value        | Description                                                                                                              |
-|-----------------------------|--------------|--------------------------------------------------------------------------------------------------------------------------|
-| **--work-package-size=**    | _1..n_       | size of work packages (defaults to _5_)                                                                                  |
-| **--ignore-result**         |              | disables store of simulation resuls                                                                                      |
-| **--dht**                   |              | enabling DHT usage (defaults to _OFF_)                                                                                   |
-| **--dht-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_)                                          |
+| **-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                    |
 
 #### Additions to `dht-snaps`
 
@@ -181,28 +210,106 @@ Following values can be set:
 ### Example: Running from scratch
 
 We will continue the above example and start a simulation with
-`dolo_diffu_inner.R`. As transport a simple fixed-coefficient diffusion is used.
-It's a 2D, 100x100 grid, simulating 10 time steps. To start the simulation with
-4 processes `cd` into your previously installed POET-dir
-`<POET_INSTALL_DIR>/bin` and run:
+*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:
 
 ```sh
-mpirun -n 4 ./poet ../share/poet/bench/dolo/dolo_diffu_inner.R/ output
+cp ../share/poet/barite/barite_het* .
+mpirun -n 4 ./poet barite_het_rt.R barite_het.qs2 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 ../share/poet/bench/dolo/dolo_diffu_inner.R/ output
+mpirun -n 4 ./poet --dht --dht-snaps=2 barite_het_rt.R barite_het.qs2 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
@@ -211,3 +318,44 @@ 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.

R_lib/CMakeLists.txt

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

R_lib/ai_surrogate_model.R

@@ -0,0 +1,75 @@
+## 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"))
+}

R_lib/init_r_lib.R

@@ -0,0 +1,112 @@
+### 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)
+}

R_lib/kin_r_library.R

@@ -1,5 +1,3 @@
-## Time-stamp: "Last modified 2023-08-15 11:58:23 delucia"
-
 ### Copyright (C) 2018-2023 Marco De Lucia, Max Luebke (GFZ Potsdam)
 ###
 ### POET is free software; you can redistribute it and/or modify it under the
@@ -15,279 +13,199 @@
 ### this program; if not, write to the Free Software Foundation, Inc., 51
 ### Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 
-
-## Simple function to check file extension. It is needed to check if
-## the GridFile is SUM (MUFITS format) or rds/RData
-FileExt <- function(x) {
-    pos <- regexpr("\\.([[:alnum:]]+)$", x)
-    ifelse(pos > -1L, substring(x, pos + 1L), "")
-}
-
-master_init <- function(setup) {
-    msgm("Process with rank 0 reading GRID properties")
-
+master_init <- function(setup, out_dir, init_field) {
     ## Setup the directory where we will store the results
-    verb <- FALSE
-    if (local_rank == 0) {
-        verb <- TRUE ## verbosity loading MUFITS results
-        if (!dir.exists(fileout)) {
-            dir.create(fileout)
-            msgm("created directory ", fileout)
-        } else {
-            msgm("dir ", fileout, " already exists, I will overwrite!")
-        }
-        if (!exists("store_result")) {
-            msgm("store_result doesn't exist!")
-        } else {
-            msgm("store_result is ", store_result)
-        }
+    if (!dir.exists(out_dir)) {
+        dir.create(out_dir)
+        msgm("created directory ", out_dir)
     } else {
-        
+        msgm("dir ", out_dir, " already exists, I will overwrite!")
+    }
+    if (is.null(setup$store_result)) {
+        msgm("store_result doesn't exist!")
+    } else {
+        msgm("store_result is ", setup$store_result)
     }
 
     setup$iter <- 1
-    setup$maxiter <- setup$iterations
     setup$timesteps <- setup$timesteps
+    setup$maxiter <- length(setup$timesteps)
+    setup$iterations <- setup$maxiter
     setup$simulation_time <- 0
 
+    dgts <- as.integer(ceiling(log10(setup$maxiter)))
+    ## string format to use in sprintf
+    fmt <- paste0("%0", dgts, "d")
+
     if (is.null(setup[["store_result"]])) {
         setup$store_result <- TRUE
     }
-    
+
     if (setup$store_result) {
+        init_field_out <- paste0(out_dir, "/iter_", sprintf(fmt = fmt, 0), ".", setup$out_ext)
+        init_field <- data.frame(init_field, check.names = FALSE)
+        SaveRObj(x = init_field, path = init_field_out)
+        msgm("Stored initial field in ", init_field_out)
         if (is.null(setup[["out_save"]])) {
             setup$out_save <- seq(1, setup$iterations)
         }
     }
-    
+
+    setup$out_dir <- out_dir
+
     return(setup)
 }
 
 ## This function, called only by master, stores on disk the last
 ## calculated time step if store_result is TRUE and increments the
 ## iteration counter
-master_iteration_end <- function(setup) {
+master_iteration_end <- function(setup, state_T, state_C) {
     iter <- setup$iter
+    # print(iter)
     ## max digits for iterations
-    dgts <-  as.integer(ceiling(log10(setup$iterations + 1)))
-    ## string format to use in sprintf 
+    dgts <- as.integer(ceiling(log10(setup$maxiter + 1)))
+    ## string format to use in sprintf
     fmt <- paste0("%0", dgts, "d")
-    
+
     ## Write on disk state_T and state_C after every iteration
     ## comprised in setup$out_save
     if (setup$store_result) {
         if (iter %in% setup$out_save) {
-            nameout <- paste0(fileout, "/iter_", sprintf(fmt=fmt, iter), ".rds")
-            info <- list(
-                tr_req_dt = as.integer(setup$req_dt)
-                ## tr_allow_dt = setup$allowed_dt,
-                ## tr_inniter = as.integer(setup$inniter)
+            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
             )
-            saveRDS(list(
-                T = setup$state_T, C = setup$state_C,
-                simtime = as.integer(setup$simtime),
-                tr_info = info
-            ), file = nameout)
+
+            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, ">")
         }
     }
-    msgm("done iteration", iter, "/", setup$maxiter)
+    ## Add last time step to simulation time
+    setup$simulation_time <- setup$simulation_time + setup$timesteps[iter]
+
+    ## msgm("done iteration", iter, "/", length(setup$timesteps))
     setup$iter <- setup$iter + 1
     return(setup)
 }
 
-## function for the workers to compute chemistry through PHREEQC
-slave_chemistry <- function(setup, data) {
-    base <- setup$base
-    first <- setup$first
-    prop <- setup$prop
-    immobile <- setup$immobile
-    kin <- setup$kin
-    ann <- setup$ann
-    
-    iter <- setup$iter
-    timesteps <- setup$timesteps
-    dt <- timesteps[iter]
-    
-    state_T <- data ## not the global field, but the work-package
-    
-    ## treat special H+/pH, e-/pe cases
-    state_T <- RedModRphree::Act2pH(state_T)
-    
-    ## reduction of the problem
-    if (setup$reduce) {
-        reduced <- ReduceStateOmit(state_T, omit = setup$ann)
-    } else {
-        reduced <- state_T
-    }
-    
-    ## form the PHREEQC input script for the current work package
-    inplist <- SplitMultiKin(
-        data = reduced, procs = 1, base = base, first = first,
-        ann = ann, prop = prop, minerals = immobile, kin = kin, dt = dt
-    )
-    
-    ## if (local_rank==1 & iter==1)
-    ##         RPhreeWriteInp("FirstInp", inplist)
-    
-    tmpC <- RunPQC(inplist, procs = 1, second = TRUE)
-    
-    ## recompose after the reduction
-    if (setup$reduce) {
-        state_C <- RecomposeState(tmpC, reduced)
-    } else {
-        state_C <- tmpC
-    }
-    
-    ## the next line is needed since we don't need all columns of
-    ## PHREEQC output
-    return(state_C[, prop])
-}
-
-## This function, called by master
-master_chemistry <- function(setup, data) {
-    state_T <- setup$state_T
-    
-    msgm(" chemistry iteration", setup$iter)
-    
-    ## treat special H+/pH, e-/pe cases
-    state_T <- RedModRphree::Act2pH(state_T)
-    
-    ## reduction of the problem
-    if (setup$reduce) {
-        reduced <- ReduceStateOmit(state_T, omit = setup$ann)
-    } else {
-        reduced <- state_T
-    }
-    
-    ## inject data from workers
-    res_C <- data
-    
-    rownames(res_C) <- NULL
-    
-    ## print(res_C)
-    
-    if (nrow(res_C) > nrow(reduced)) {
-        res_C <- res_C[seq(2, nrow(res_C), by = 2), ]
-    }
-    
-    ## recompose after the reduction
-    if (setup$reduce) {
-        state_C <- RecomposeState(res_C, reduced)
-    } else {
-        state_C <- res_C
-    }
-    
-    setup$state_C <- state_C
-    setup$reduced <- reduced
-    
-    return(setup)
-}
-
-
-## Adapted version for "reduction"
-ReduceStateOmit <- function(data, omit = NULL, sign = 6) {
-    require(mgcv)
-    
-    rem <- colnames(data)
-    if (is.list(omit)) {
-        indomi <- match(names(omit), colnames(data))
-        datao <- data[, -indomi]
-    } else {
-        datao <- data
-    }
-    
-    datao <- signif(datao, sign)
-    red <- mgcv::uniquecombs(datao)
-    inds <- attr(red, "index")
-    now <- ncol(red)
-    
-    
-    ## reattach the omitted column(s)
-    ## FIXME: control if more than one ann is present
-    if (is.list(omit)) {
-        red <- cbind(red, rep(data[1, indomi], nrow(red)))
-
-        colnames(red)[now + 1] <- names(omit)
-
-        ret <- red[, colnames(data)]
-    } else {
-        ret <- red
-    }
-    rownames(ret) <- NULL
-    attr(ret, "index") <- inds
-    return(ret)
-}
-
-
 
 ## Attach the name of the calling function to the message displayed on
 ## R's stdout
 msgm <- function(...) {
-    if (local_rank == 0) {
-        fname <- as.list(sys.call(-1))[[1]]
-        prefix <- paste0("R: ", fname, " ::")
-        cat(paste(prefix, ..., "\n"))
-    }
+    prefix <- paste0("R: ")
+    cat(paste(prefix, ..., "\n"))
     invisible()
 }
 
 
-## Function called by master R process to store on disk all relevant
-## parameters for the simulation
-StoreSetup <- function(setup) {
-
-    to_store <- vector(mode = "list", length = 4)
-    ## names(to_store) <- c("Sim", "Flow", "Transport", "Chemistry", "DHT")
-    names(to_store) <- c("Sim", "Transport", "DHT", "Cmdline")
-    
-    ## read the setup R file, which is sourced in kin.cpp
-    tmpbuff <- file(filesim, "r")
-    setupfile <- readLines(tmpbuff)
-    close.connection(tmpbuff)
-    
-    to_store$Sim <- setupfile
-    
-    ## to_store$Flow <- list(
-    ##     snapshots  = setup$snapshots,
-    ##     gridfile   = setup$gridfile,
-    ##     phase      = setup$phase,
-    ##     density    = setup$density,
-    ##     dt_differ  = setup$dt_differ,
-    ##     prolong    = setup$prolong,
-    ##     maxiter    = setup$maxiter,
-    ##     saved_iter = setup$iter_output,
-    ##     out_save   = setup$out_save )
-
-    to_store$Transport <- setup$diffusion
-
-    ## to_store$Chemistry <- list(
-    ##    nprocs   = n_procs,
-    ##    wp_size  = work_package_size,
-    ##    base     = setup$base,
-    ##    first    = setup$first,
-    ##    init     = setup$initsim,
-    ##    db       = db,
-    ##    kin      = setup$kin,
-    ##    ann      = setup$ann)
-
-    if (dht_enabled) {
-        to_store$DHT <- list(
-            enabled   = dht_enabled,
-            log       = dht_log
-            ## signif    = dht_final_signif,
-            ## proptype  = dht_final_proptype
-        )
-    } else {
-        to_store$DHT <- FALSE
-    }
-
-    to_store$Cmdline <- commandArgs(trailingOnly=FALSE)
-    saveRDS(to_store, file = paste0(fileout, "/setup.rds"))
-    msgm("initialization stored in ", paste0(fileout, "/setup.rds"))
-}
-
 GetWorkPackageSizesVector <- function(n_packages, package_size, len) {
-    ids <- rep(1:n_packages, times=package_size, each = 1)[1:len]
+    ids <- rep(1:n_packages, times = package_size, each = 1)[1:len]
     return(as.integer(table(ids)))
 }
+
+
+## Handler to read R objs from binary files using either builtin
+## readRDS(), qs::qread() or qs2::qs_read() based on file extension
+ReadRObj <- function(path) {
+    ## code borrowed from tools::file_ext()
+    pos <- regexpr("\\.([[:alnum:]]+)$", path)
+    extension <- ifelse(pos > -1L, substring(path, pos + 1L), "")
+
+    switch(extension,
+        rds = readRDS(path),
+        qs  = qs::qread(path),
+        qs2 = qs2::qs_read(path)
+    )
+}
+
+## Handler to store R objs to binary files using either builtin
+## saveRDS() or qs::qsave() based on file extension
+SaveRObj <- function(x, path) {
+    ## msgm("Storing to", path)
+    ## code borrowed from tools::file_ext()
+    pos <- regexpr("\\.([[:alnum:]]+)$", path)
+    extension <- ifelse(pos > -1L, substring(path, pos + 1L), "")
+
+    switch(extension,
+        rds = saveRDS(object = x, file = path),
+        qs  = qs::qsave(x = x, file = path),
+        qs2 = qs2::qs_save(object = x, file = path)
+    )
+}
+
+
+######## Old relic code
+
+## ## Function called by master R process to store on disk all relevant
+## ## parameters for the simulation
+## StoreSetup <- function(setup, filesim, out_dir) {
+##     to_store <- vector(mode = "list", length = 4)
+##     ## names(to_store) <- c("Sim", "Flow", "Transport", "Chemistry", "DHT")
+##     names(to_store) <- c("Sim", "Transport", "DHT", "Cmdline")
+
+##     ## read the setup R file, which is sourced in kin.cpp
+##     tmpbuff <- file(filesim, "r")
+##     setupfile <- readLines(tmpbuff)
+##     close.connection(tmpbuff)
+
+##     to_store$Sim <- setupfile
+
+##     ## to_store$Flow <- list(
+##     ##     snapshots  = setup$snapshots,
+##     ##     gridfile   = setup$gridfile,
+##     ##     phase      = setup$phase,
+##     ##     density    = setup$density,
+##     ##     dt_differ  = setup$dt_differ,
+##     ##     prolong    = setup$prolong,
+##     ##     maxiter    = setup$maxiter,
+##     ##     saved_iter = setup$iter_output,
+##     ##     out_save   = setup$out_save )
+
+##     to_store$Transport <- setup$diffusion
+
+##     ## to_store$Chemistry <- list(
+##     ##    nprocs   = n_procs,
+##     ##    wp_size  = work_package_size,
+##     ##    base     = setup$base,
+##     ##    first    = setup$first,
+##     ##    init     = setup$initsim,
+##     ##    db       = db,
+##     ##    kin      = setup$kin,
+##     ##    ann      = setup$ann)
+
+##     if (dht_enabled) {
+##         to_store$DHT <- list(
+##             enabled   = dht_enabled,
+##             log       = dht_log
+##             ## signif    = dht_final_signif,
+##             ## proptype  = dht_final_proptype
+##         )
+##     } else {
+##         to_store$DHT <- FALSE
+##     }
+
+##     if (dht_enabled) {
+##         to_store$DHT <- list(
+##             enabled   = dht_enabled,
+##             log       = dht_log
+##             # signif    = dht_final_signif,
+##             # proptype  = dht_final_proptype
+##         )
+##     } else {
+##         to_store$DHT <- FALSE
+##     }
+
+##     saveRDS(to_store, file = paste0(fileout, "/setup.rds"))
+##     msgm("initialization stored in ", paste0(fileout, "/setup.rds"))
+## }

app/CMakeLists.txt

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

app/poet.cpp

@@ -1,369 +0,0 @@
-/*
-** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
-** Potsdam)
-**
-** Copyright (C) 2018-2022 Marco De Lucia, Max Luebke (GFZ Potsdam)
-**
-** POET is free software; you can redistribute it and/or modify it under the
-** terms of the GNU General Public License as published by the Free Software
-** Foundation; either version 2 of the License, or (at your option) any later
-** version.
-**
-** POET is distributed in the hope that it will be useful, but WITHOUT ANY
-** WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
-** A PARTICULAR PURPOSE. See the GNU General Public License for more details.
-**
-** You should have received a copy of the GNU General Public License along with
-** this program; if not, write to the Free Software Foundation, Inc., 51
-** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
-*/
-
-#include <Rcpp.h>
-#include <cstdint>
-#include <cstdlib>
-#include <poet/ChemistryModule.hpp>
-#include <poet/DiffusionModule.hpp>
-#include <poet/Grid.hpp>
-#include <poet/Macros.hpp>
-#include <poet/RInsidePOET.hpp>
-#include <poet/SimParams.hpp>
-
-#include <cstring>
-#include <iostream>
-#include <string>
-#include <vector>
-
-#include <mpi.h>
-#include <poet.h>
-
-using namespace std;
-using namespace poet;
-using namespace Rcpp;
-
-poet::ChemistryModule::SingleCMap DFToHashMap(const Rcpp::DataFrame &df) {
-  std::unordered_map<std::string, double> out_map;
-  vector<string> col_names = Rcpp::as<vector<string>>(df.names());
-
-  for (const auto &name : col_names) {
-    double val = df[name.c_str()];
-    out_map.insert({name, val});
-  }
-
-  return out_map;
-}
-
-// HACK: this is a step back as the order and also the count of fields is
-// predefined, but it will change in the future
-void writeFieldsToR(RInside &R, const Field &trans, const Field &chem) {
-  R["TMP"] = Rcpp::wrap(trans.AsVector());
-  R["TMP_PROPS"] = Rcpp::wrap(trans.GetProps());
-  R.parseEval(std::string(
-      "mysetup$state_T <- setNames(data.frame(matrix(TMP, nrow=" +
-      std::to_string(trans.GetRequestedVecSize()) + ")), TMP_PROPS)"));
-
-  R["TMP"] = Rcpp::wrap(chem.AsVector());
-  R["TMP_PROPS"] = Rcpp::wrap(chem.GetProps());
-  R.parseEval(std::string(
-      "mysetup$state_C <- setNames(data.frame(matrix(TMP, nrow=" +
-      std::to_string(chem.GetRequestedVecSize()) + ")), TMP_PROPS)"));
-}
-
-void set_chem_parameters(poet::ChemistryModule &chem, uint32_t wp_size,
-                         const std::string &database_path) {
-  chem.SetErrorHandlerMode(1);
-  chem.SetComponentH2O(false);
-  chem.SetRebalanceFraction(0.5);
-  chem.SetRebalanceByCell(true);
-  chem.UseSolutionDensityVolume(false);
-  chem.SetPartitionUZSolids(false);
-
-  // Set concentration units
-  // 1, mg/L; 2, mol/L; 3, kg/kgs
-  chem.SetUnitsSolution(2);
-  // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
-  chem.SetUnitsPPassemblage(1);
-  // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
-  chem.SetUnitsExchange(1);
-  // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
-  chem.SetUnitsSurface(1);
-  // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
-  chem.SetUnitsGasPhase(1);
-  // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
-  chem.SetUnitsSSassemblage(1);
-  // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
-  chem.SetUnitsKinetics(1);
-
-  // Set representative volume
-  std::vector<double> rv;
-  rv.resize(wp_size, 1.0);
-  chem.SetRepresentativeVolume(rv);
-
-  // Set initial porosity
-  std::vector<double> por;
-  por.resize(wp_size, 1);
-  chem.SetPorosity(por);
-
-  // Set initial saturation
-  std::vector<double> sat;
-  sat.resize(wp_size, 1.0);
-  chem.SetSaturation(sat);
-
-  // Load database
-  chem.LoadDatabase(database_path);
-}
-
-inline double RunMasterLoop(SimParams &params, RInside &R,
-                            const GridParams &g_params, uint32_t nxyz_master) {
-
-  DiffusionParams d_params{R};
-  DiffusionModule diffusion(d_params, g_params);
-  /* Iteration Count is dynamic, retrieving value from R (is only needed by
-   * master for the following loop) */
-  uint32_t maxiter = R.parseEval("mysetup$iterations");
-
-  double sim_time = .0;
-
-  ChemistryModule chem(nxyz_master, params.getNumParams().wp_size, maxiter,
-                       params.getChemParams(), MPI_COMM_WORLD);
-
-  set_chem_parameters(chem, nxyz_master, params.getChemParams().database_path);
-  chem.RunInitFile(params.getChemParams().input_script);
-
-  poet::ChemistryModule::SingleCMap init_df = DFToHashMap(d_params.initial_t);
-  chem.initializeField(diffusion.getField());
-
-  if (params.getNumParams().print_progressbar) {
-    chem.setProgressBarPrintout(true);
-  }
-
-  /* SIMULATION LOOP */
-
-  double dSimTime{0};
-  for (uint32_t iter = 1; iter < maxiter + 1; iter++) {
-    double start_t = MPI_Wtime();
-    uint32_t tick = 0;
-    // cout << "CPP: Evaluating next time step" << endl;
-    // R.parseEvalQ("mysetup <- master_iteration_setup(mysetup)");
-
-    double dt = Rcpp::as<double>(
-        R.parseEval("mysetup$timesteps[" + std::to_string(iter) + "]"));
-
-    //  cout << "CPP: Next time step is " << dt << "[s]" << endl;
-    MSG("Next time step is " + std::to_string(dt) + " [s]");
-
-    /* displaying iteration number, with C++ and R iterator */
-    MSG("Going through iteration " + std::to_string(iter));
-    MSG("R's $iter: " +
-        std::to_string((uint32_t)(R.parseEval("mysetup$iter"))) +
-        ". Iteration");
-
-    /* run transport */
-    // TODO: transport to diffusion
-    diffusion.simulate(dt);
-
-    chem.getField().update(diffusion.getField());
-
-    MSG("Chemistry step");
-
-    chem.SetTimeStep(dt);
-    chem.RunCells();
-
-    writeFieldsToR(R, diffusion.getField(), chem.GetField());
-    diffusion.getField().update(chem.GetField());
-
-    R["req_dt"] = dt;
-    R["simtime"] = (sim_time += dt);
-
-    R.parseEval("mysetup$req_dt <- req_dt");
-    R.parseEval("mysetup$simtime <- simtime");
-
-    // MDL master_iteration_end just writes on disk state_T and
-    // state_C after every iteration if the cmdline option
-    // --ignore-results is not given (and thus the R variable
-    // store_result is TRUE)
-    R.parseEvalQ("mysetup <- master_iteration_end(setup=mysetup)");
-
-    MSG("End of *coupling* iteration " + std::to_string(iter) + "/" +
-        std::to_string(maxiter));
-    MSG();
-
-    // MPI_Barrier(MPI_COMM_WORLD);
-    double end_t = MPI_Wtime();
-    dSimTime += end_t - start_t;
-  } // END SIMULATION LOOP
-
-  R.parseEvalQ("profiling <- list()");
-
-  R["simtime_chemistry"] = chem.GetChemistryTime();
-  R.parseEvalQ("profiling$simtime_chemistry <- simtime_chemistry");
-
-  R["chemistry_loop"] = chem.GetMasterLoopTime();
-  R.parseEvalQ("profiling$chemistry_loop <- chemistry_loop");
-
-  R["chemistry_sequential"] = chem.GetMasterSequentialTime();
-  R.parseEvalQ("profiling$simtime_sequential <- chemistry_sequential");
-
-  R["idle_master"] = chem.GetMasterIdleTime();
-  R.parseEvalQ("profiling$idle_master <- idle_master");
-
-  R["idle_worker"] = Rcpp::wrap(chem.GetWorkerIdleTimings());
-  R.parseEvalQ("profiling$idle_worker <- idle_worker");
-
-  R["phreeqc_time"] = Rcpp::wrap(chem.GetWorkerPhreeqcTimings());
-  R.parseEvalQ("profiling$phreeqc <- phreeqc_time");
-
-  R["simtime_transport"] = diffusion.getTransportTime();
-  R.parseEvalQ("profiling$simtime_transport <- simtime_transport");
-
-  // R["phreeqc_count"] = phreeqc_counts;
-  // R.parseEvalQ("profiling$phreeqc_count <- phreeqc_count");
-
-  if (params.getChemParams().use_dht) {
-    R["dht_hits"] = Rcpp::wrap(chem.GetWorkerDHTHits());
-    R.parseEvalQ("profiling$dht_hits <- dht_hits");
-    R["dht_evictions"] = Rcpp::wrap(chem.GetWorkerDHTEvictions());
-    R.parseEvalQ("profiling$dht_evictions <- dht_evictions");
-    R["dht_get_time"] = Rcpp::wrap(chem.GetWorkerDHTGetTimings());
-    R.parseEvalQ("profiling$dht_get_time <- dht_get_time");
-    R["dht_fill_time"] = Rcpp::wrap(chem.GetWorkerDHTFillTimings());
-    R.parseEvalQ("profiling$dht_fill_time <- dht_fill_time");
-  }
-  if (params.getChemParams().use_interp) {
-    R["interp_w"] = Rcpp::wrap(chem.GetWorkerInterpolationWriteTimings());
-    R.parseEvalQ("profiling$interp_write <- interp_w");
-    R["interp_r"] = Rcpp::wrap(chem.GetWorkerInterpolationReadTimings());
-    R.parseEvalQ("profiling$interp_read <- interp_r");
-    R["interp_g"] = Rcpp::wrap(chem.GetWorkerInterpolationGatherTimings());
-    R.parseEvalQ("profiling$interp_gather <- interp_g");
-    R["interp_fc"] =
-        Rcpp::wrap(chem.GetWorkerInterpolationFunctionCallTimings());
-    R.parseEvalQ("profiling$interp_function_calls <- interp_fc");
-    R["interp_calls"] = Rcpp::wrap(chem.GetWorkerInterpolationCalls());
-    R.parseEvalQ("profiling$interp_calls <- interp_calls");
-    R["interp_cached"] = Rcpp::wrap(chem.GetWorkerPHTCacheHits());
-    R.parseEvalQ("profiling$interp_cached <- interp_cached");
-  }
-
-  chem.MasterLoopBreak();
-  diffusion.end();
-
-  return dSimTime;
-}
-
-int main(int argc, char *argv[]) {
-  double dSimTime, sim_end;
-
-  int world_size, world_rank;
-
-  MPI_Init(&argc, &argv);
-
-  MPI_Comm_size(MPI_COMM_WORLD, &world_size);
-
-  MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);
-
-  RInsidePOET &R = RInsidePOET::getInstance();
-
-  if (world_rank == 0) {
-    MSG("Running POET version " + std::string(poet_version));
-  }
-
-  if (world_rank > 0) {
-    {
-      SimParams params(world_rank, world_size);
-      int pret = params.parseFromCmdl(argv, R);
-
-      if (pret == poet::PARSER_ERROR) {
-        MPI_Finalize();
-        return EXIT_FAILURE;
-      } else if (pret == poet::PARSER_HELP) {
-        MPI_Finalize();
-        return EXIT_SUCCESS;
-      }
-
-      // ChemistryModule worker(nxyz, nxyz, MPI_COMM_WORLD);
-      ChemistryModule worker = poet::ChemistryModule::createWorker(
-          MPI_COMM_WORLD, params.getChemParams());
-      set_chem_parameters(worker, worker.GetWPSize(),
-                          params.getChemParams().database_path);
-
-      worker.WorkerLoop();
-    }
-
-    MPI_Barrier(MPI_COMM_WORLD);
-
-    MSG("finished, cleanup of process " + std::to_string(world_rank));
-
-    MPI_Finalize();
-
-    return EXIT_SUCCESS;
-  }
-
-  /*Loading Dependencies*/
-  // TODO: kann raus
-  std::string r_load_dependencies = "source('../R_lib/kin_r_library.R');";
-  R.parseEvalQ(r_load_dependencies);
-
-  SimParams params(world_rank, world_size);
-  int pret = params.parseFromCmdl(argv, R);
-
-  if (pret == poet::PARSER_ERROR) {
-    MPI_Finalize();
-    return EXIT_FAILURE;
-  } else if (pret == poet::PARSER_HELP) {
-    MPI_Finalize();
-    return EXIT_SUCCESS;
-  }
-
-  MSG("RInside initialized on process " + std::to_string(world_rank));
-
-  R.parseEvalQ("mysetup <- setup");
-  // if (world_rank == 0) { // get timestep vector from
-  // grid_init function ... //
-  std::string master_init_code = "mysetup <- master_init(setup=setup)";
-  R.parseEval(master_init_code);
-
-  GridParams g_params(R);
-
-  params.initVectorParams(R);
-
-  // MDL: store all parameters
-  if (world_rank == 0) {
-    MSG("Calling R Function to store calling parameters");
-    R.parseEvalQ("StoreSetup(setup=mysetup)");
-  }
-
-  if (world_rank == 0) {
-    MSG("Init done on process with rank " + std::to_string(world_rank));
-  }
-
-  // MPI_Barrier(MPI_COMM_WORLD);
-
-  uint32_t nxyz_master = (world_size == 1 ? g_params.total_n : 1);
-
-  dSimTime = RunMasterLoop(params, R, g_params, nxyz_master);
-
-  MSG("finished simulation loop");
-
-  MSG("start timing profiling");
-
-  R["simtime"] = dSimTime;
-  R.parseEvalQ("profiling$simtime <- simtime");
-
-  string r_vis_code;
-  r_vis_code = "saveRDS(profiling, file=paste0(fileout,'/timings.rds'));";
-  R.parseEval(r_vis_code);
-
-  MSG("Done! Results are stored as R objects into <" + params.getOutDir() +
-      "/timings.rds>");
-
-  MPI_Barrier(MPI_COMM_WORLD);
-
-  MSG("finished, cleanup of process " + std::to_string(world_rank));
-  MPI_Finalize();
-
-  if (world_rank == 0) {
-    MSG("done, bye!");
-  }
-
-  exit(0);
-}

app/poet.h.in

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

apps/CMakeLists.txt

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

apps/initializer/main.cpp

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

bench/CMakeLists.txt

@@ -1,3 +1,43 @@
+function(ADD_BENCH_TARGET TARGET POET_BENCH_LIST RT_FILES OUT_PATH)
+    set(bench_install_dir share/poet/${OUT_PATH})
+
+    # create empty list 
+    set(OUT_FILES_LIST "")
+
+    foreach(BENCH_FILE ${${POET_BENCH_LIST}})
+        get_filename_component(BENCH_NAME ${BENCH_FILE} NAME_WE)
+        set(OUT_FILE ${CMAKE_CURRENT_BINARY_DIR}/${BENCH_NAME})
+        set(OUT_FILE_EXT ${OUT_FILE}.qs2)
+
+        add_custom_command(
+            OUTPUT ${OUT_FILE_EXT}
+            COMMAND $<TARGET_FILE:poet_init> -n ${OUT_FILE} -s ${CMAKE_CURRENT_SOURCE_DIR}/${BENCH_FILE}
+            COMMENT "Running poet_init on ${BENCH_FILE}"
+            DEPENDS poet_init ${CMAKE_CURRENT_SOURCE_DIR}/${BENCH_FILE}
+            VERBATIM
+            COMMAND_EXPAND_LISTS
+        )
+
+        list(APPEND OUT_FILES_LIST ${OUT_FILE_EXT})
+
+    endforeach(BENCH_FILE ${${POET_BENCH_LIST}})
+        
+    add_custom_target(
+        ${TARGET}
+        DEPENDS ${OUT_FILES_LIST})
+            
+    install(FILES ${OUT_FILES_LIST} DESTINATION ${bench_install_dir})
+
+    # install all ADD_FILES to the same location
+    install(FILES ${${RT_FILES}} DESTINATION ${bench_install_dir})
+endfunction()
+
+
+# define target name 
+set(BENCHTARGET benchmarks)
+
+add_custom_target(${BENCHTARGET} ALL)
+
+add_subdirectory(barite)
 add_subdirectory(dolo)
 add_subdirectory(surfex)
-add_subdirectory(barite)

bench/barite/CMakeLists.txt

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

bench/barite/README.org

@@ -18,9 +18,19 @@ mpirun -np 4 ./poet --interp  barite_interp_eval.R barite_results
 
 * List of Files
 
-- =barite.R=: POET input script for a 20x20 simulation grid
-- =barite_interp_eval.R=: POET input script for a 400x200 simulation
+- =barite_het.R=: POET input script with homogeneous zones for a 5x2 simulation
+  grid 
+- =barite_200.R=: POET input script for a 200x200 simulation
   grid
+- =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

bench/barite/barite.R

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

bench/barite/barite.pqi

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

bench/barite/barite_200.R

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

bench/barite/barite_200_rt.R

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

bench/barite/barite_200ai_surrogate_input_script.R

@@ -0,0 +1,48 @@
+## 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)
+}

bench/barite/barite_50ai.R

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

bench/barite/barite_50ai_rt.R

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

bench/barite/barite_50ai_surr_mdl.R

@@ -0,0 +1,90 @@
+## 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)
+}

bench/barite/barite_het.R

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

bench/barite/barite_het.pqi

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

bench/barite/barite_het_rt.R

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

bench/barite/barite_interp_eval.R

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

bench/dolo/CMakeLists.txt

@@ -1,9 +1,18 @@
-install(FILES
-    dolo_diffu_inner.R
-    dolo_diffu_inner_large.R
-    dolo_inner.pqi
-    dolo_interp_long.R
-    phreeqc_kin.dat
-DESTINATION
-    share/poet/bench/dolo
+set(bench_files
+    dolo_inner_large.R
+    dolo_interp.R
 )
+
+set(runtime_files
+    dolo_inner_large_rt.R
+    dolo_interp_rt.R
+)
+
+ADD_BENCH_TARGET(
+    dolo_bench
+    bench_files 
+    runtime_files
+    "dolo"
+)
+
+add_dependencies(${BENCHTARGET} dolo_bench)

bench/dolo/Eval.R

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

bench/dolo/README.org

@@ -18,16 +18,13 @@ mpirun -np 4 ./poet --dht --interp dolo_interp_long.R dolo_interp_long_res
 
 * List of Files
 
-- =dolo_diffu_inner.R=: POET input script for a 100x100 simulation
-  grid
-- =dolo_interp_long.R=: POET input script for a 400x200 simulation
+- =dolo_interp.R=: POET input script for a 400x200 simulation
   grid
+- =dolo_diffu_inner_large.R=: POET input script for a 400x200
+  simulation grid
 - =phreeqc_kin.dat=: PHREEQC database containing the kinetic expressions
   for dolomite and celestite, stripped down from =phreeqc.dat=
 - =dol.pqi=: PHREEQC input script for the chemical system
-# - =dolo.R=: POET input script for a 20x20 simulation grid
-# - =dolo_diffu_inner_large.R=: POET input script for a 400x200
-#   simulation grid
 
 * Chemical system
 

bench/dolo/dol.pqi

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

bench/dolo/dolo_diffu_inner.R

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

bench/dolo/dolo_diffu_inner_large.R

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

bench/dolo/dolo_inner.pqi

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

bench/dolo/dolo_inner_large.R

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

bench/dolo/dolo_inner_large_rt.R

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

bench/dolo/dolo_interp.R

@@ -0,0 +1,131 @@
+rows <- 400
+cols <- 200
+
+grid_def <- matrix(2, nrow = rows, ncol = cols)
+
+# Define grid configuration for POET model
+grid_setup <- list(
+    pqc_in_file = "./dol.pqi",
+    pqc_db_file = "./phreeqc_kin.dat", # Path to the database file for Phreeqc
+    grid_def = grid_def, # Definition of the grid, containing IDs according to the Phreeqc input script
+    grid_size = c(2.5, 5), # Size of the grid in meters
+    constant_cells = c() # IDs of cells with constant concentration
+)
+
+bound_def_we <- list(
+    "type" = rep("constant", rows),
+    "sol_id" = rep(1, rows),
+    "cell" = seq(1, rows)
+)
+
+bound_def_ns <- list(
+    "type" = rep("constant", cols),
+    "sol_id" = rep(1, cols),
+    "cell" = seq(1, cols)
+)
+
+diffusion_setup <- list(
+    boundaries = list(
+        "W" = bound_def_we,
+        "E" = bound_def_we,
+        "N" = bound_def_ns,
+        "S" = bound_def_ns
+    ),
+    inner_boundaries = list(
+        "row" = floor(rows / 2),
+        "col" = floor(cols / 2),
+        "sol_id" = c(3)
+    ),
+    alpha_x = 1e-6,
+    alpha_y = 1e-6
+)
+
+check_sign_cal_dol_dht <- function(old, new) {
+    # if ((old["Calcite"] == 0) != (new["Calcite"] == 0)) {
+    #     return(TRUE)
+    # }
+    # if ((old["Dolomite"] == 0) != (new["Dolomite"] == 0)) {
+    #     return(TRUE)
+    # }
+    return(FALSE)
+}
+
+# fuzz_input_dht_keys <- function(input) {
+#     dht_species <- c(
+#         "H" = 3,
+#         "O" = 3,
+#         "Charge" = 3,
+#         "C" = 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
+)

bench/dolo/dolo_interp_long.R

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

bench/dolo/dolo_interp_rt.R

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

bench/fgcs/20241211_README.tex

@@ -0,0 +1,102 @@
+% 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:

bench/fgcs/EvalFGCS.R

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

bench/fgcs/README.org

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

bench/fgcs/barite_fgcs_2.R

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

bench/fgcs/barite_fgcs_2.pqi

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

bench/fgcs/barite_fgcs_2_rt.R

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

bench/surfex/CMakeLists.txt

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

bench/surfex/ExBase.pqi

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

bench/surfex/PoetEGU_surfex_500.R

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

bench/surfex/PoetEGU_surfex_500_rt.R

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

bench/surfex/README.org

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

bench/surfex/SurfExBase.pqi

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

bench/surfex/SurfexEGU.pqi

@@ -0,0 +1,108 @@
+##  Time-stamp: "Last modified 2024-04-12 10:59:59 delucia"
+## KNOBS
+##  -logfile            false
+##  -iterations         10000
+##  -convergence_tolerance  1E-12
+##  -step_size          2
+##  -pe_step_size       2
+
+SOLUTION 1                     ## Porewater composition Opalinus Clay, WITHOUT radionuclides, AFTER EQUI_PHASES
+  pe           -2.627          ## Eh = -227 mV, Value from Bossart & Thury (2008)-> PC borehole measurement 2003, Eh still decreasing  
+  density      1.01583         ## kg/dm³ = g/cm³
+  temp         13              ## mean temperature Mont Terri, Bossart & Thury (2008), calculations performed for 25°C
+  units        mol/kgw
+  ## Mean composition   
+  pH           7.064                    
+  Na           2.763e-01
+  Cl           3.228e-01  charge 
+  S(6)         1.653e-02  as SO4
+  Ca           2.173e-02 
+  Mg           1.740e-02
+  K            1.902e-03 
+  Sr           4.520e-04 
+  Fe           1.435e-04                 
+  U            2.247e-09                
+
+SURFACE 1 Opalinus Clay, clay minerals 
+  ## calculated with rho_b=2.2903 kg/dm³, poro=0.1662
+  ## 1 dm³ = 13.565641 kg_sed/kg_pw
+  -equil 1                ## equilibrate with solution 1 
+  -sites_units density    ## set unit for binding site density to sites/nm2
+  -donnan 4.9e-10         ## calculated after Wigger & Van Loon (2018) for ionic strength after equilibration with minerales for pCO2=2.2 log10 bar    
+
+  # surface  density SSA (m2/g)  mass (g/kgw) 
+  Kln_aOH  1.155   11.         3798.4     ## Kaolinite 28 wt% (aluminol and silanol sites) 
+  Kln_siOH 1.155               
+  Ill_sOH  0.05    100.        4205.35    ## Illite 31 wt% (weak und strong binding sites) 
+  Ill_wOH  2.26                           ## 2 % strong binding sites  
+  Mll_sOH  0.05    100.        813.94     ## Montmorillonite = smektite = 6 wt% (weak und strong binding sites)            
+  Mll_wOH  2.26                           ## 2 % strong binding sites
+
+EXCHANGE 1 Exchanger, only illite+montmorillonite 
+  ## Illite = 0.225 eq/kg_rock, Montmorillonit = 0.87 eq/kg_rock 
+ -equil 1         ## equilibrate with solution 1
+  Z     0.9462    ## = Illite
+  Y     0.70813   ## = Montmorillonite
+END
+
+SOLUTION 2                     ##  Porewater composition Opalinus Clay, WITHOUT radionuclides, AFTER EQUI_PHASES
+  pe           -2.627          ## Eh = -227 mV, Value from Bossart & Thury (2008)-> PC borehole measurement 2003, Eh still decreasing  
+  density      1.01583         ## kg/dm³ = g/cm³
+  temp         13              ## mean temperature Mont Terri, Bossart & Thury (2008), calculations performed for 25°C
+  units        mol/kgw
+  ## Mean composition   
+  pH           7.064                    
+  Na           2.763e-01
+  Cl           3.228e-01  charge 
+  S(6)         1.653e-02  as SO4
+  Ca           2.173e-02 
+  Mg           1.740e-02
+  K            1.902e-03 
+  Sr           4.520e-04 
+  Fe           1.435e-04                 
+  U            2.247e-09                
+
+SURFACE 2 Opalinus Clay, clay minerals 
+  -equil 2                ## equilibrate with solution 2 
+  -sites_units density    ## set unit for binding site density to
+			  ## sites/nm2
+  -donnan 4.9e-10         ## calculated after Wigger & Van Loon (2018)
+			  ## for ionic strength after equilibration
+			  ## with minerales for pCO2=2.2 log10 bar
+
+  ## surface  density SSA (m2/g)  mass (g/kgw) 
+  Kln_aOH  1.155   11.         2798.4     ## Kaolinite 28 wt% (aluminol and silanol sites) 
+  Kln_siOH 1.155               
+  Ill_sOH  0.05    100.        1205.35    ## Illite 31 wt% (weak und strong binding sites) 
+  Ill_wOH  2.26                           ## 2 % strong binding sites  
+  Mll_sOH  0.05    100.        113.94     ## Montmorillonite = smektite = 6 wt% (weak und strong binding sites)            
+  Mll_wOH  2.26                           ## 2 % strong binding sites
+
+EXCHANGE 2 Exchanger, only illite+montmorillonite 
+  ## Illite = 0.225 eq/kg_rock, Montmorillonit = 0.87 eq/kg_rock 
+ -equil 2         ## equilibrate with solution 1
+  Z     0.5       ## = Illite
+  Y     0.2       ## = Montmorillonite
+END
+
+SOLUTION 3
+  pe           -2.627          ## Eh = -227 mV, Value from Bossart & Thury (2008)-> PC borehole measurement 2003, Eh still decreasing  
+  density      1.01583         ## kg/dm³ = g/cm³
+  temp         13              ## mean temperature Mont Terri, Bossart & Thury (2008), calculations performed for 25°C
+  units        mol/kgw
+  ## Mean composition   
+  pH           7.064                    
+  Na           3.763e-01
+  Cl           4.228e-01  charge 
+  S(6)         1.653e-02  as SO4
+  Ca           2.173e-02 
+  Mg           1.740e-02
+  K            1.902e-03 
+  Sr           4.520e-04 
+  Fe           1.435e-04                 
+  U            1e-6            
+  C            1.991e-03 
+END
+
+RUN_CELLS
+END

bench/surfex/ex.R

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

bench/surfex/ex_rt.R

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

bench/surfex/surfex.R

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

bench/surfex/surfex_rt.R

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

docs/CMakeLists.txt

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

docs/Input_Scripts.md

@@ -1,86 +0,0 @@
-# Input Scripts
-
-In the following the expected schemes of the input scripts is described.
-Therefore, each section of the input script gets its own chapter. All sections
-should return a `list` as results, which are concatenated to one setup list at
-the end of the file. All values must have the same name in order to get parsed
-by POET.
-
-## Grid initialization
-
-| name           | type           | description                                                           |
-|----------------|----------------|-----------------------------------------------------------------------|
-| `n_cells`      | Numeric Vector | Number of cells in each direction                                     |
-| `s_cells`      | Numeric Vector | Spatial resolution of grid in each direction                          |
-| `type`         | String         | Type of initialization, can be set to *scratch*, *phreeqc* or *rds*   |
-
-## Diffusion parameters
-
-| name           | type                 | description                               |
-|----------------|----------------------|-------------------------------------------|
-| `init`         | Named Numeric Vector | Initial state for each diffused species   |
-| `vecinj`       | Data Frame           | Defining all boundary conditions row wise |
-| `vecinj_inner` | List of Triples      | Inner boundaries                          |
-| `vecinj_index` | List of 4 elements   | Ghost nodes boundary conditions           |
-| `alpha`        | Named Numeric Vector | Constant alpha for each species           |
-
-### Remark on boundary conditions
-
-Each boundary condition should be defined in `vecinj` as a data frame, where one
-row holds one boundary condition.
-
-To define inner (constant) boundary conditions, use a list of triples in
-`vecinj_inner`, where each triples is defined by $(i,x,y)$. $i$ is defining the
-boundary condition, referencing to the row in `vecinj`. $x$ and $y$ coordinates
-then defining the position inside the grid. 
-
-Ghost nodes are set by `vecinj_index` which is a list containing boundaries for
-each celestial direction (**important**: named by `N, E, S, W`). Each direction
-is a numeric vector, also representing a row index of the `vecinj` data frame
-for each ghost node, starting at the left-most and upper cell respectively. By
-setting the boundary condition to $0$, the ghost node is set as closed boundary.
-
-#### Example
-
-Suppose you have a `vecinj` data frame defining 2 boundary conditions and a grid
-consisting of $10 \times 10$ grid cells. Grid cell $(1,1)$ should be set to the
-first boundary condition and $(5,6)$ to the second. Also, all boundary
-conditions for the ghost nodes should be closed. Except the southern boundary,
-which should be set to the first boundary condition injection. The following
-setup describes how to setup your initial script, where `n` and `m` are the
-grids cell count for each direction ($n = m = 10$):
-
-```R
-vecinj_inner <- list (
-  l1 = c(1, 1, 1),
-  l2 = c(2, 5, 6)
-)
-
-vecinj_index <- list(
-  "N" = rep(0, n),
-  "E" = rep(0, m),
-  "S" = rep(1, n),
-  "W" = rep(0, m)
-)
-```
-
-## Chemistry parameters
-
-| name           | type         | description                                                                      |
-|----------------|--------------|----------------------------------------------------------------------------------|
-| `database`     | String       | Path to the Phreeqc database                                                     |
-| `input_script` | String       | Path the the Phreeqc input script                                                |
-| `dht_species`  | Named Vector | Indicates significant digits to use for each species for DHT rounding.           |
-| `pht_species`  | Named Vector | Indicates significant digits to use for each species for Interpolation rounding. |
-
-## Final setup
-
-| name           | type           | description                                                |
-|----------------|----------------|------------------------------------------------------------|
-| `grid`         | List           | Grid parameter list                                        |
-| `diffusion`    | List           | Diffusion parameter list                                   |
-| `chemistry`    | List           | Chemistry parameter list                                   |
-| `iterations`   | Numeric Value  | Count of iterations                                        |
-| `timesteps`    | Numeric Vector | $\Delta t$ to use for specific iteration                   |
-| `store_result` | Boolean        | Indicates if results should be stored                      |
-| `out_save`     | Numeric Vector | *optional:* At which iteration the states should be stored |

docs/Output.md

@@ -35,34 +35,50 @@ 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 |
-| simtime\_transport | measured time in *transport* subroutine                                    |
-| simtime\_chemistry | measured time in *chemistry* subroutine (actual parallelized part)         |
+| Value     | Description                                                                |
+| --------- | -------------------------------------------------------------------------- |
+| simtime   | time spent in whole simulation loop without any initialization and cleanup |
+| chemistry | measured time in *chemistry* subroutine                                    |
+| diffusion | measured time in *diffusion* subroutine                                    |
 
-### chemistry subsetting
+### Chemistry subsetting
 
-If running parallel there are also measured timings which are subsets of
-*simtime\_​chemistry*.
+| Value         | Description                                               |
+| ------------- | --------------------------------------------------------- |
+| simtime       | overall runtime of chemistry                              |
+| loop          | time spent in send/recv loop of master                    |
+| sequential    | sequential part of the master (e.g. shuffling field)      |
+| idle\_master  | idling time of the master waiting for workers             |
+| idle\_worker  | idling time (waiting for work from master) of the workers |
+| phreeqc\_time | accumulated times for Phreeqc calls of every worker       |
 
-| Value                 | Description                                               |
-|-----------------------|-----------------------------------------------------------|
-| chemistry\_loop       | time spent in send/​recv loop of master                    |
-| chemistry\_sequential | sequential part of master chemistry                       |
-| idle\_master          | idling time (waiting for any free worker) of the master   |
-| idle\_worker          | idling time (waiting for work from master) of the workers |
-| phreeqc\_time         | accumulated times for Phreeqc calls of every worker       |
-
-### DHT usage {#DHT-usage}
+#### DHT usage
 
 If running in parallel and with activated DHT, two more timings and also
 some profiling about the DHT usage are given:
 
 | Value           | Description                                             |
-|-----------------|---------------------------------------------------------|
-| dht\_fill\_time | time to write data to DHT                               |
-| dht\_get\_time  | time to retreive data from DHT                          |
-| dh\_hits        | count of data points retrieved from DHT                 |
-| dht\_miss       | count of misses/count of data points written to DHT     |
+| --------------- | ------------------------------------------------------- |
+| dht\_hits       | count of data points retrieved from DHT                 |
 | dht\_evictions  | count of data points evicted by another write operation |
+| dht\_get\_time  | time to retreive data from DHT                          |
+| dht\_fill\_time | time to write data to DHT                               |
+
+#### Interpolation
+
+If using interpolation, the following values are given:
+
+| Value          | Description                                                           |
+| -------------- | --------------------------------------------------------------------- |
+| interp\_w      | time spent to write to PHT                                            |
+| interp\_r      | time spent to read from DHT/PHT/Cache                                 |
+| interp\_g      | time spent to gather results from DHT                                 |
+| interp\_fc     | accumulated time spent in interpolation function call                 |
+| interp\_calls  | count of interpolations                                               |
+| interp\_cached | count of interpolation data sets, which where cached in the local map |
+
+### Diffusion subsetting
+
+| Value     | Description                                |
+| --------- | ------------------------------------------ |
+| simtime   | overall runtime of diffusion               |

ext/doctest

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

ext/litephreeqc

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

ext/phreeqcrm

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

ext/tug

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

include/poet/Grid.hpp

@@ -1,116 +0,0 @@
-/*
-** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
-** Potsdam)
-**
-** Copyright (C) 2018-2023 Marco De Lucia, Max Luebke (GFZ Potsdam)
-**
-** POET is free software; you can redistribute it and/or modify it under the
-** terms of the GNU General Public License as published by the Free Software
-** Foundation; either version 2 of the License, or (at your option) any later
-** version.
-**
-** POET is distributed in the hope that it will be useful, but WITHOUT ANY
-** WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
-** A PARTICULAR PURPOSE. See the GNU General Public License for more details.
-**
-** You should have received a copy of the GNU General Public License along with
-** this program; if not, write to the Free Software Foundation, Inc., 51
-** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
-*/
-
-#ifndef GRID_H
-#define GRID_H
-
-#include "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

include/poet/RInsidePOET.hpp

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

include/poet/SimParams.hpp

@@ -1,273 +0,0 @@
-/*
-** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
-** Potsdam)
-**
-** Copyright (C) 2018-2022 Marco De Lucia, Max Luebke (GFZ Potsdam)
-**
-** POET is free software; you can redistribute it and/or modify it under the
-** terms of the GNU General Public License as published by the Free Software
-** Foundation; either version 2 of the License, or (at your option) any later
-** version.
-**
-** POET is distributed in the hope that it will be useful, but WITHOUT ANY
-** WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
-** A PARTICULAR PURPOSE. See the GNU General Public License for more details.
-**
-** You should have received a copy of the GNU General Public License along with
-** this program; if not, write to the Free Software Foundation, Inc., 51
-** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
-*/
-
-#ifndef PARSER_H
-#define PARSER_H
-
-#include <cstdint>
-#include <optional>
-#include <string>
-#include <string_view>
-#include <unordered_map>
-#include <vector>
-
-#include "DataStructures.hpp"
-#include "Macros.hpp"
-#include "RInsidePOET.hpp"
-#include "argh.hpp" // Argument handler https://github.com/adishavit/argh
-#include <RInside.h>
-#include <Rcpp.h>
-// BSD-licenced
-
-/** Standard DHT Size. Defaults to 1 GB (1000 MB) */
-constexpr uint32_t DHT_SIZE_PER_PROCESS_MB = 1.5E3;
-/** Standard work package size */
-#define WORK_PACKAGE_SIZE_DEFAULT 32
-
-namespace poet {
-
-enum { PARSER_OK, PARSER_ERROR, PARSER_HELP };
-
-/**
- * @brief Defining all simulation parameters
- *
- */
-typedef struct {
-  /** Count of processes in MPI_COMM_WORLD */
-  int world_size;
-  /** rank of proces in MPI_COMM_WORLD */
-  int world_rank;
-  /** indicates if DHT should be used */
-  bool dht_enabled;
-  /** apply logarithm to key before rounding */
-  bool dht_log;
-  /** indicates if timestep dt differs between iterations */
-  bool dt_differ;
-  /** Indicates, when a DHT snapshot should be written */
-  int dht_snaps;
-  /** <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;
-  /** indicating whether the progress bar during chemistry simulation should be
-   * printed or not */
-  bool print_progressbar;
-
-  bool interp_enabled;
-} t_simparams;
-
-using GridParams = struct s_GridParams {
-  std::array<uint32_t, 2> n_cells;
-  std::array<double, 2> s_cells;
-  std::uint8_t dim;
-  std::uint32_t total_n;
-
-  std::string type;
-
-  Rcpp::DataFrame init_df;
-  std::string input_script;
-  std::string database_path;
-
-  std::vector<std::string> props;
-
-  s_GridParams(RInside &R);
-};
-
-using DiffusionParams = struct s_DiffusionParams {
-  Rcpp::DataFrame initial_t;
-
-  Rcpp::NumericVector alpha;
-  Rcpp::List vecinj_inner;
-
-  Rcpp::DataFrame vecinj;
-  Rcpp::DataFrame vecinj_index;
-
-  s_DiffusionParams(RInside &R);
-};
-
-struct ChemistryParams {
-  std::string database_path;
-  std::string input_script;
-
-  bool use_dht;
-  std::uint64_t dht_size;
-  int dht_snaps;
-  std::string dht_file;
-  std::string dht_outdir;
-  NamedVector<std::uint32_t> dht_signifs;
-
-  bool use_interp;
-  std::uint64_t pht_size;
-  std::uint32_t pht_max_entries;
-  std::uint32_t interp_min_entries;
-  NamedVector<std::uint32_t> pht_signifs;
-
-  struct Chem_Hook_Functions {
-    RHookFunction<bool> dht_fill;
-    RHookFunction<std::vector<double>> dht_fuzz;
-    RHookFunction<std::vector<std::size_t>> interp_pre;
-    RHookFunction<bool> interp_post;
-  } hooks;
-
-  void initFromR(RInsidePOET &R);
-};
-
-/**
- * @brief Reads information from program arguments and R runtime
- *
- * Providing functions to initialize parameters of the simulation using command
- * line parameters and parameters from the R runtime. This class will also parse
- * arguments from the commandline and decides if argument is known or unknown.
- *
- * Stores and distribute current simulation parameters at any time.
- *
- */
-class SimParams {
-public:
-  /**
-   * @brief Construct a new SimParams object
-   *
-   * With all given parameters a new instance of this class will be created.
-   *
-   * @param world_rank Rank of process inside MPI_COMM_WORLD
-   * @param world_size Size of communicator MPI_COMM_WORLD
-   */
-  SimParams(int world_rank, int world_size);
-
-  /**
-   * @brief Parse program arguments
-   *
-   * This is done by the argh.h library.
-   *
-   * First, the function will check if there is a flag 'help' or 'h'. If this is
-   * the case a help message is printed and the function will return with
-   * PARSER_HELP.
-   *
-   * Second, if there are not 2 positional arguments an error will be printed to
-   * stderr and the function returns with PARSER_ERROR.
-   *
-   * Then all given program parameters and flags will be read and checked, if
-   * there are known by validateOptions. A list of all unknown options might be
-   * returned, printed out and the function will return with PARSER_ERROR.
-   * Oterhwise the function continuos.
-   *
-   * Now all program arguments will be stored inside t_simparams struct, printed
-   * out and the function returns with PARSER_OK.
-   *
-   * Also, all parsed agruments are distributed to the R runtime.
-   *
-   * @param argv Argument value of the program
-   * @param R Instantiated R runtime
-   * @return int Returns with 0 if no error occured, otherwise value less than 0
-   * is returned.
-   */
-  int parseFromCmdl(char *argv[], RInsidePOET &R);
-
-  /**
-   * @brief Init std::vector values
-   *
-   * This will initialize dht_signif_vector and dht_prop_type_vector internally
-   * depending on whether vectors are defined by R-Simulation file or not.
-   * 'init_chemistry' must be run beforehand.
-   *
-   * @param R R runtime
-   */
-  void initVectorParams(RInside &R);
-
-  /**
-   * @brief Get the numerical params struct
-   *
-   * Returns a struct which contains all numerical or boolean simulation
-   * parameters.
-   *
-   * @return t_simparams Parameter struct
-   */
-  auto getNumParams() const { return this->simparams; };
-
-  /**
-   * @brief Get the DHT_Signif_Vector
-   *
-   * Returns a vector indicating which significant values are used for each
-   * variable of a grid cell.
-   *
-   * @return std::vector<int> Vector of integers containing information about
-   * significant digits
-   */
-  auto getDHTSignifVector() const { return this->dht_signif_vector; };
-
-  auto getPHTSignifVector() const { return this->pht_signif_vector; };
-
-  auto getPHTBucketSize() const { return this->pht_bucket_size; };
-  auto getPHTMinEntriesNeeded() const { return this->pht_min_entries_needed; };
-
-  /**
-   * @brief Get the filesim name
-   *
-   * Returns a string containing the absolute path to a R file defining the
-   * simulation.
-   *
-   * @return std::string Absolute path to R file
-   */
-  auto getFilesim() const { return this->filesim; };
-
-  /**
-   * @brief Get the output directory
-   *
-   * Returns the name of an absolute path where all output files should be
-   * stored.
-   *
-   * @return std::string Absolute path to output path
-   */
-  auto getOutDir() const { return this->out_dir; };
-
-  const auto &getChemParams() const { return chem_params; }
-
-private:
-  std::list<std::string> validateOptions(argh::parser cmdl);
-
-  const std::set<std::string> flaglist{"ignore-result", "dht", "P", "progress",
-                                       "interp"};
-  const std::set<std::string> paramlist{
-      "work-package-size", "dht-strategy",
-      "dht-size",          "dht-snaps",
-      "dht-file",          "interp-size",
-      "interp-min",        "interp-bucket-entries"};
-
-  t_simparams simparams;
-
-  std::vector<uint32_t> dht_signif_vector;
-  std::vector<uint32_t> pht_signif_vector;
-
-  uint32_t pht_bucket_size;
-  uint32_t pht_min_entries_needed;
-
-  std::string filesim;
-  std::string out_dir;
-
-  ChemistryParams chem_params;
-};
-} // namespace poet
-#endif // PARSER_H

include/poet/argh.hpp

@@ -1,459 +0,0 @@
-/*
-** 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));
-   }
-}
-
-

include/poet/enums.hpp

@@ -1,10 +0,0 @@
-#ifndef ENUMS_H_
-#define ENUMS_H_
-
-namespace poet {
-enum DHT_PROP_TYPES { DHT_TYPE_DEFAULT, DHT_TYPE_CHARGE, DHT_TYPE_TOTAL };
-
-enum CHEMISTRY_OUT_SOURCE { CHEM_PQC, CHEM_DHT, CHEM_INTERP };
-} // namespace poet
-
-#endif // ENUMS_H_

src/Base/Macros.hpp

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

src/Base/RInsidePOET.hpp

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

src/CMakeLists.txt

@@ -1,25 +1,101 @@
-file(GLOB_RECURSE poet_lib_SRC
-  CONFIGURE_DEPENDS
-  "*.cpp" "*.c")
+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_library(poet_lib ${poet_lib_SRC})
-target_include_directories(poet_lib PUBLIC ${PROJECT_SOURCE_DIR}/include)
-target_link_libraries(poet_lib PUBLIC
-  MPI::MPI_CXX ${MATH_LIBRARY} RRuntime PhreeqcRM tug)
+set(POET_TUG_APPROACH "Implicit" CACHE STRING "tug numerical approach to use")
+set_property(CACHE POET_TUG_APPROACH PROPERTY STRINGS "Implicit" "Explicit")
 
-target_compile_definitions(poet_lib 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(poet_lib PRIVATE DHT_STATISTICS)
+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()
 
-option(POET_PHT_ADDITIONAL_INFO "Enables additional information in the PHT" OFF)
+target_include_directories(POETLib PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}")
+target_link_libraries(
+  POETLib 
+  PUBLIC RRuntime
+  PUBLIC IPhreeqcPOET
+  PUBLIC tug
+  PUBLIC MPI::MPI_C
+)
 
-if (POET_PHT_ADDITIONAL_INFO)
-  target_compile_definitions(poet_lib PRIVATE POET_PHT_ADD)
-endif()
+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)

src/Chemistry/ChemistryDefs.hpp

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

src/Chemistry/ChemistryModule.cpp

@@ -0,0 +1,324 @@
+#include "ChemistryModule.hpp"
+
+#include "PhreeqcEngine.hpp"
+#include "PhreeqcMatrix.hpp"
+#include "PhreeqcRunner.hpp"
+#include "SurrogateModels/DHT_Wrapper.hpp"
+#include "SurrogateModels/Interpolation.hpp"
+
+#include <algorithm>
+#include <cassert>
+#include <cmath>
+#include <cstdint>
+#include <memory>
+#include <mpi.h>
+#include <string>
+#include <vector>
+
+std::vector<double>
+inverseDistanceWeighting(const std::vector<std::int32_t> &to_calc,
+                         const std::vector<double> &from,
+                         const std::vector<std::vector<double>> &input,
+                         const std::vector<std::vector<double>> &output) {
+  std::vector<double> results = from;
+
+  // const std::uint32_t buffer_size = input.size() + 1;
+  // double buffer[buffer_size];
+  // double from_rescaled;
+
+  const std::uint32_t data_set_n = input.size();
+  double rescaled[to_calc.size()][data_set_n + 1];
+  double weights[data_set_n];
+
+  // rescaling over all key elements
+  for (std::uint32_t key_comp_i = 0; key_comp_i < to_calc.size();
+       key_comp_i++) {
+    const auto output_comp_i = to_calc[key_comp_i];
+
+    // rescale input between 0 and 1
+    for (std::uint32_t point_i = 0; point_i < data_set_n; point_i++) {
+      rescaled[key_comp_i][point_i] = input[point_i][key_comp_i];
+    }
+
+    rescaled[key_comp_i][data_set_n] = from[output_comp_i];
+
+    const double min = *std::min_element(rescaled[key_comp_i],
+                                         rescaled[key_comp_i] + data_set_n + 1);
+    const double max = *std::max_element(rescaled[key_comp_i],
+                                         rescaled[key_comp_i] + data_set_n + 1);
+
+    for (std::uint32_t point_i = 0; point_i < data_set_n; point_i++) {
+      rescaled[key_comp_i][point_i] =
+          ((max - min) != 0
+               ? (rescaled[key_comp_i][point_i] - min) / (max - min)
+               : 0);
+    }
+    rescaled[key_comp_i][data_set_n] =
+        ((max - min) != 0 ? (from[output_comp_i] - min) / (max - min) : 0);
+  }
+
+  // calculate distances for each data set
+  double inv_sum = 0;
+  for (std::uint32_t point_i = 0; point_i < data_set_n; point_i++) {
+    double distance = 0;
+    for (std::uint32_t key_comp_i = 0; key_comp_i < to_calc.size();
+         key_comp_i++) {
+      distance += std::pow(
+          rescaled[key_comp_i][point_i] - rescaled[key_comp_i][data_set_n], 2);
+    }
+
+    weights[point_i] = distance != 0 ? 1 / std::sqrt(distance) : 0;
+    assert(!std::isnan(weights[point_i]));
+    inv_sum += weights[point_i];
+  }
+
+  assert(!std::isnan(inv_sum));
+
+  // actual interpolation
+  // bool has_h = false;
+  // bool has_o = false;
+
+  for (std::uint32_t key_comp_i = 0; key_comp_i < to_calc.size();
+       key_comp_i++) {
+    const auto output_comp_i = to_calc[key_comp_i];
+    double key_delta = 0;
+
+    // if (interp_i == 0) {
+    //   has_h = true;
+    // }
+
+    // if (interp_i == 1) {
+    //   has_o = true;
+    // }
+
+    for (std::uint32_t j = 0; j < data_set_n; j++) {
+      key_delta += weights[j] * output[j][output_comp_i];
+    }
+
+    key_delta /= inv_sum;
+
+    results[output_comp_i] = from[output_comp_i] + key_delta;
+    assert(!std::isnan(results[output_comp_i]));
+  }
+
+  return results;
+}
+
+poet::ChemistryModule::ChemistryModule(
+    uint32_t wp_size_, const InitialList::ChemistryInit chem_params,
+    MPI_Comm communicator)
+    : group_comm(communicator), wp_size(wp_size_), params(chem_params) {
+  MPI_Comm_rank(communicator, &comm_rank);
+  MPI_Comm_size(communicator, &comm_size);
+
+  this->is_sequential = comm_size == 1;
+  this->is_master = comm_rank == 0;
+
+  this->n_cells = chem_params.total_grid_cells;
+
+  if (!is_master) {
+    PhreeqcMatrix pqc_mat =
+        PhreeqcMatrix(chem_params.database, chem_params.pqc_script,
+                      chem_params.with_h0_o0, chem_params.with_redox);
+
+    this->pqc_runner =
+        std::make_unique<PhreeqcRunner>(pqc_mat.subset(chem_params.pqc_ids));
+  }
+}
+
+poet::ChemistryModule::~ChemistryModule() {
+  if (dht) {
+    delete dht;
+  }
+}
+
+void poet::ChemistryModule::initializeDHT(
+    uint32_t size_mb, const NamedVector<std::uint32_t> &key_species,
+    bool has_het_ids) {
+  constexpr uint32_t MB_FACTOR = 1E6;
+
+  MPI_Comm dht_comm;
+
+  if (this->is_master) {
+    MPI_Comm_split(this->group_comm, MPI_UNDEFINED, this->comm_rank, &dht_comm);
+    return;
+  }
+
+  if (!this->is_master) {
+
+    MPI_Comm_split(this->group_comm, 1, this->comm_rank, &dht_comm);
+
+    auto map_copy = key_species;
+
+    if (key_species.empty()) {
+      std::vector<std::uint32_t> default_signif(
+          this->prop_count, DHT_Wrapper::DHT_KEY_SIGNIF_DEFAULT);
+      map_copy = NamedVector<std::uint32_t>(this->prop_names, default_signif);
+    }
+
+    auto key_indices = parseDHTSpeciesVec(key_species, this->prop_names);
+
+    if (this->dht) {
+      delete this->dht;
+    }
+
+    const std::uint64_t dht_size = size_mb * MB_FACTOR;
+
+    this->dht = new DHT_Wrapper(dht_comm, dht_size, map_copy, key_indices,
+                                this->prop_names, params.hooks,
+                                this->prop_count, interp_enabled, has_het_ids);
+    this->dht->setBaseTotals(base_totals.at(0), base_totals.at(1));
+  }
+}
+
+inline std::vector<std::int32_t> poet::ChemistryModule::parseDHTSpeciesVec(
+    const NamedVector<std::uint32_t> &key_species,
+    const std::vector<std::string> &to_compare) const {
+  std::vector<int32_t> species_indices;
+  species_indices.reserve(key_species.size());
+
+  const auto test = key_species.getNames();
+
+  for (const auto &species : key_species.getNames()) {
+    auto it = std::find(to_compare.begin(), to_compare.end(), species);
+    if (it == to_compare.end()) {
+      species_indices.push_back(DHT_Wrapper::DHT_KEY_INPUT_CUSTOM);
+      continue;
+    }
+    const std::uint32_t index = it - to_compare.begin();
+    species_indices.push_back(index);
+  }
+
+  return species_indices;
+}
+
+void poet::ChemistryModule::BCastStringVec(std::vector<std::string> &io) {
+
+  if (this->is_master) {
+    int vec_size = io.size();
+    ChemBCast(&vec_size, 1, MPI_INT);
+
+    for (const auto &value : io) {
+      int buf_size = value.size() + 1;
+      ChemBCast(&buf_size, 1, MPI_INT);
+      ChemBCast(const_cast<char *>(value.c_str()), buf_size, MPI_CHAR);
+    }
+  } else {
+    int vec_size;
+    ChemBCast(&vec_size, 1, MPI_INT);
+
+    io.resize(vec_size);
+
+    for (int i = 0; i < vec_size; i++) {
+      int buf_size;
+      ChemBCast(&buf_size, 1, MPI_INT);
+      char buf[buf_size];
+      ChemBCast(buf, buf_size, MPI_CHAR);
+      io[i] = std::string{buf};
+    }
+  }
+}
+
+void poet::ChemistryModule::setDHTSnapshots(int type,
+                                            const std::string &out_dir) {
+  if (this->is_master) {
+    return;
+  }
+
+  this->dht_file_out_dir = out_dir;
+  this->dht_snaps_type = type;
+}
+
+void poet::ChemistryModule::setDHTReadFile(const std::string &input_file) {
+  if (this->is_master) {
+    return;
+  }
+
+  if (!input_file.empty()) {
+    WorkerReadDHTDump(input_file);
+  }
+}
+
+void poet::ChemistryModule::initializeInterp(
+    std::uint32_t bucket_size, std::uint32_t size_mb, std::uint32_t min_entries,
+    const NamedVector<std::uint32_t> &key_species) {
+
+  if (!this->is_master) {
+
+    constexpr uint32_t MB_FACTOR = 1E6;
+
+    assert(this->dht);
+
+    this->interp_enabled = true;
+
+    auto map_copy = key_species;
+
+    if (key_species.empty()) {
+      map_copy = this->dht->getKeySpecies();
+      for (auto i = 0; i < map_copy.size(); i++) {
+        const std::uint32_t signif =
+            static_cast<std::uint32_t>(map_copy[i]) -
+            (map_copy[i] > InterpolationModule::COARSE_DIFF
+                 ? InterpolationModule::COARSE_DIFF
+                 : 0);
+        map_copy[i] = signif;
+      }
+    }
+
+    auto key_indices =
+        parseDHTSpeciesVec(map_copy, dht->getKeySpecies().getNames());
+
+    if (this->interp) {
+      this->interp.reset();
+    }
+
+    const uint64_t pht_size = size_mb * MB_FACTOR;
+
+    interp = std::make_unique<poet::InterpolationModule>(
+        bucket_size, pht_size, min_entries, *(this->dht), map_copy, key_indices,
+        this->prop_names, this->params.hooks);
+
+    interp->setInterpolationFunction(inverseDistanceWeighting);
+  }
+}
+
+std::vector<double>
+poet::ChemistryModule::shuffleField(const std::vector<double> &in_field,
+                                    uint32_t size_per_prop, uint32_t prop_count,
+                                    uint32_t wp_count) {
+  std::vector<double> out_buffer(in_field.size());
+  uint32_t write_i = 0;
+  for (uint32_t i = 0; i < wp_count; i++) {
+    for (uint32_t j = i; j < size_per_prop; j += wp_count) {
+      for (uint32_t k = 0; k < prop_count; k++) {
+        out_buffer[(write_i * prop_count) + k] =
+            in_field[(k * size_per_prop) + j];
+      }
+      write_i++;
+    }
+  }
+  return out_buffer;
+}
+
+void poet::ChemistryModule::unshuffleField(const std::vector<double> &in_buffer,
+                                           uint32_t size_per_prop,
+                                           uint32_t prop_count,
+                                           uint32_t wp_count,
+                                           std::vector<double> &out_field) {
+  uint32_t read_i = 0;
+
+  for (uint32_t i = 0; i < wp_count; i++) {
+    for (uint32_t j = i; j < size_per_prop; j += wp_count) {
+      for (uint32_t k = 0; k < prop_count; k++) {
+        out_field[(k * size_per_prop) + j] =
+            in_buffer[(read_i * prop_count) + k];
+      }
+      read_i++;
+    }
+  }
+}
+
+void poet::ChemistryModule::set_ai_surrogate_validity_vector(
+    std::vector<int> r_vector) {
+  this->ai_surrogate_validity_vector = r_vector;
+}

src/Chemistry/ChemistryModule.hpp

@@ -1,17 +1,19 @@
-//  Time-stamp: "Last modified 2023-08-15 14:36:28 mluebke"
 
 #ifndef CHEMISTRYMODULE_H_
 #define CHEMISTRYMODULE_H_
 
-#include "DHT_Wrapper.hpp"
-#include "DataStructures.hpp"
-#include "Interpolation.hpp"
-#include "IrmResult.h"
-#include "PhreeqcRM.h"
-#include "SimParams.hpp"
+#include "DataStructures/Field.hpp"
+#include "DataStructures/NamedVector.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>
@@ -24,7 +26,7 @@ namespace poet {
  * \brief Wrapper around PhreeqcRM to provide POET specific parallelization with
  * easy access.
  */
-class ChemistryModule : public PhreeqcRM {
+class ChemistryModule {
 public:
   /**
    * Creates a new instance of Chemistry module with given grid cell count, work
@@ -41,61 +43,80 @@ public:
    * \param wp_size Count of grid cells to fill each work package at maximum.
    * \param communicator MPI communicator to distribute work in.
    */
-  ChemistryModule(uint32_t nxyz, uint32_t wp_size, std::uint32_t maxiter,
-                  const ChemistryParams &chem_param, MPI_Comm communicator);
+  ChemistryModule(uint32_t wp_size,
+                  const InitialList::ChemistryInit chem_params,
+                  MPI_Comm communicator);
 
   /**
    * Deconstructor, which frees DHT data structure if used.
    */
   ~ChemistryModule();
 
-  /**
-   * Parses the input script and extract information needed during runtime.
-   *
-   * **Only run by master**.
-   *
-   * Database must be loaded beforehand.
-   *
-   * \param input_script_path Path to input script to parse.
-   */
-  void RunInitFile(const std::string &input_script_path);
-
+  void masterSetField(Field field);
   /**
    * Run the chemical simulation with parameters set.
    */
-  void RunCells();
+  void simulate(double dt);
 
-  /**
-   * Returns the chemical field.
-   */
-  auto GetField() const { return this->chem_field; }
   /**
    * Returns all known species names, including not only aqueous species, but
    * also equilibrium, exchange, surface and kinetic reactants.
    */
-  auto GetPropNames() const { return this->prop_names; }
+  // auto GetPropNames() const { return this->prop_names; }
 
   /**
    * Return the accumulated runtime in seconds for chemical simulation.
    */
   auto GetChemistryTime() const { return this->chem_t; }
 
-  /**
-   * Create a new worker instance inside given MPI communicator.
-   *
-   * Wraps communication needed before instanciation can take place.
-   *
-   * \param communicator MPI communicator to distribute work in.
-   *
-   * \returns A worker instance with fixed work package size.
-   */
-  static ChemistryModule createWorker(MPI_Comm communicator,
-                                      const ChemistryParams &chem_param);
+  void setFilePadding(std::uint32_t maxiter) {
+    this->file_pad =
+        static_cast<std::uint8_t>(std::ceil(std::log10(maxiter + 1)));
+  }
 
-  /**
-   * Default work package size.
-   */
-  static constexpr uint32_t CHEM_DEFAULT_WP_SIZE = 5;
+  struct SurrogateSetup {
+    std::vector<std::string> prop_names;
+    std::array<double, 2> base_totals;
+    bool has_het_ids;
+
+    bool dht_enabled;
+    std::uint32_t dht_size_mb;
+    int dht_snaps;
+    std::string dht_out_dir;
+
+    bool interp_enabled;
+    std::uint32_t interp_bucket_size;
+    std::uint32_t interp_size_mb;
+    std::uint32_t interp_min_entries;
+    bool ai_surrogate_enabled;
+  };
+
+  void masterEnableSurrogates(const SurrogateSetup &setup) {
+    // FIXME: This is a hack to get the prop_names and prop_count from the setup
+    this->prop_names = setup.prop_names;
+    this->prop_count = setup.prop_names.size();
+
+    this->dht_enabled = setup.dht_enabled;
+    this->interp_enabled = setup.interp_enabled;
+    this->ai_surrogate_enabled = setup.ai_surrogate_enabled;
+
+    this->base_totals = setup.base_totals;
+
+    if (this->dht_enabled || this->interp_enabled) {
+      this->initializeDHT(setup.dht_size_mb, this->params.dht_species,
+                          setup.has_het_ids);
+
+      if (setup.dht_snaps != DHT_SNAPS_DISABLED) {
+        this->setDHTSnapshots(setup.dht_snaps, setup.dht_out_dir);
+      }
+    }
+
+    if (this->interp_enabled) {
+      this->initializeInterp(setup.interp_bucket_size, setup.interp_size_mb,
+                             setup.interp_min_entries,
+                             this->params.interp_species);
+    }
+  }
 
   /**
    * Intended to alias input parameters for grid initialization with a single
@@ -118,17 +139,6 @@ public:
     DHT_SNAPS_ITEREND       //!< output snapshots after each iteration
   };
 
-  /**
-   * **This function has to be run!**
-   *
-   * Merge initial values from existing module with the chemistry module and set
-   * according internal variables.
-   *
-   * \param other Field to merge chemistry with. Most likely it is something
-   * like the diffusion field.
-   */
-  void initializeField(const Field &other);
-
   /**
    * **Only called by workers!** Start the worker listening loop.
    */
@@ -223,6 +233,11 @@ public:
     this->print_progessbar = enabled;
   };
 
+  /**
+   *  **Master only** Set the ai surrogate validity vector from R
+   */
+  void set_ai_surrogate_validity_vector(std::vector<int> r_vector);
+
   std::vector<uint32_t> GetWorkerInterpolationCalls() const;
 
   std::vector<double> GetWorkerInterpolationWriteTimings() const;
@@ -232,9 +247,12 @@ public:
 
   std::vector<uint32_t> GetWorkerPHTCacheHits() const;
 
+  std::vector<int> ai_surrogate_validity_vector;
+
 protected:
   void initializeDHT(uint32_t size_mb,
-                     const NamedVector<std::uint32_t> &key_species);
+                     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);
 
@@ -243,12 +261,9 @@ protected:
                         const NamedVector<std::uint32_t> &key_species);
 
   enum {
-    CHEM_INIT,
-    CHEM_WP_SIZE,
-    CHEM_INIT_SPECIES,
+    CHEM_FIELD_INIT,
     CHEM_DHT_ENABLE,
     CHEM_DHT_SIGNIF_VEC,
-    CHEM_DHT_PROP_TYPE_VEC,
     CHEM_DHT_SNAPS,
     CHEM_DHT_READ_FILE,
     CHEM_IP_ENABLE,
@@ -256,7 +271,8 @@ protected:
     CHEM_IP_SIGNIF_VEC,
     CHEM_WORK_LOOP,
     CHEM_PERF,
-    CHEM_BREAK_MAIN_LOOP
+    CHEM_BREAK_MAIN_LOOP,
+    CHEM_AI_BCAST_VALIDITY
   };
 
   enum { LOOP_WORK, LOOP_END };
@@ -295,7 +311,7 @@ protected:
   using worker_list_t = std::vector<struct worker_info_s>;
   using workpointer_t = std::vector<double>::iterator;
 
-  void MasterRunParallel();
+  void MasterRunParallel(double dt);
   void MasterRunSequential();
 
   void MasterSendPkgs(worker_list_t &w_list, workpointer_t &work_pointer,
@@ -321,8 +337,8 @@ protected:
   void WorkerPerfToMaster(int type, const struct worker_s &timings);
   void WorkerMetricsToMaster(int type);
 
-  IRM_RESULT WorkerRunWorkPackage(WorkPackage &work_package, double dSimTime,
-                                  double dTimestep);
+  void WorkerRunWorkPackage(WorkPackage &work_package, double dSimTime,
+                            double dTimestep);
 
   std::vector<uint32_t> CalculateWPSizesVector(uint32_t n_cells,
                                                uint32_t wp_size) const;
@@ -345,7 +361,7 @@ protected:
   bool is_sequential;
   bool is_master;
 
-  uint32_t wp_size{CHEM_DEFAULT_WP_SIZE};
+  uint32_t wp_size;
   bool dht_enabled{false};
   int dht_snaps_type{DHT_SNAPS_DISABLED};
   std::string dht_file_out_dir;
@@ -355,6 +371,8 @@ protected:
   bool interp_enabled{false};
   std::unique_ptr<poet::InterpolationModule> interp;
 
+  bool ai_surrogate_enabled{false};
+
   static constexpr uint32_t BUFFER_OFFSET = 5;
 
   inline void ChemBCast(void *buf, int count, MPI_Datatype datatype) const {
@@ -373,7 +391,7 @@ protected:
 
   bool print_progessbar{false};
 
-  std::uint32_t file_pad;
+  std::uint8_t file_pad{1};
 
   double chem_t{0.};
 
@@ -383,11 +401,9 @@ protected:
 
   Field chem_field;
 
-  static constexpr int MODULE_COUNT = 5;
+  const InitialList::ChemistryInit params;
 
-  const ChemistryParams &params;
-
-  std::array<std::uint32_t, MODULE_COUNT> speciesPerModule{};
+  std::unique_ptr<PhreeqcRunner> pqc_runner;
 };
 } // namespace poet
 

src/Chemistry/MasterFunctions.cpp

@@ -1,11 +1,9 @@
-#include "PhreeqcRM.h"
-#include "poet/ChemistryModule.hpp"
+#include "ChemistryModule.hpp"
 
 #include <algorithm>
 #include <cstddef>
 #include <cstdint>
 #include <mpi.h>
-#include <stdexcept>
 #include <vector>
 
 std::vector<uint32_t>
@@ -161,6 +159,20 @@ std::vector<uint32_t> poet::ChemistryModule::GetWorkerPHTCacheHits() const {
   return ret;
 }
 
+inline std::vector<int> shuffleVector(const std::vector<int> &in_vector,
+                                      uint32_t size_per_prop,
+                                      uint32_t wp_count) {
+  std::vector<int> out_buffer(in_vector.size());
+  uint32_t write_i = 0;
+  for (uint32_t i = 0; i < wp_count; i++) {
+    for (uint32_t j = i; j < size_per_prop; j += wp_count) {
+      out_buffer[write_i] = in_vector[j];
+      write_i++;
+    }
+  }
+  return out_buffer;
+}
+
 inline std::vector<double> shuffleField(const std::vector<double> &in_field,
                                         uint32_t size_per_prop,
                                         uint32_t prop_count,
@@ -249,8 +261,10 @@ inline void poet::ChemistryModule::MasterSendPkgs(
       send_buffer[end_of_wp + 2] = dt;
       // current time of simulation (age) in seconds
       send_buffer[end_of_wp + 3] = this->simtime;
-      // placeholder for work_package_count
-      send_buffer[end_of_wp + 4] = 0.;
+      // current work package start location in field
+      uint32_t wp_start_index = std::accumulate(wp_sizes_vector.begin(), std::next(wp_sizes_vector.begin(), count_pkgs), 0);
+      send_buffer[end_of_wp + 4] = wp_start_index;
+
 
       /* ATTENTION Worker p has rank p+1 */
       // MPI_Send(send_buffer, end_of_wp + BUFFER_OFFSET, MPI_DOUBLE, p + 1,
@@ -307,45 +321,46 @@ inline void poet::ChemistryModule::MasterRecvPkgs(worker_list_t &w_list,
   }
 }
 
-void poet::ChemistryModule::RunCells() {
+void poet::ChemistryModule::simulate(double dt) {
   double start_t{MPI_Wtime()};
   if (this->is_sequential) {
     MasterRunSequential();
     return;
   }
 
-  MasterRunParallel();
+  MasterRunParallel(dt);
   double end_t{MPI_Wtime()};
   this->chem_t += end_t - start_t;
 }
 
 void poet::ChemistryModule::MasterRunSequential() {
-  std::vector<double> shuffled_field =
-      shuffleField(chem_field.AsVector(), n_cells, prop_count, 1);
+  // std::vector<double> shuffled_field =
+  //     shuffleField(chem_field.AsVector(), n_cells, prop_count, 1);
 
-  std::vector<std::vector<double>> input;
-  for (std::size_t i = 0; i < n_cells; i++) {
-    input.push_back(
-        std::vector<double>(shuffled_field.begin() + (i * prop_count),
-                            shuffled_field.begin() + ((i + 1) * prop_count)));
-  }
+  // std::vector<std::vector<double>> input;
+  // for (std::size_t i = 0; i < n_cells; i++) {
+  //   input.push_back(
+  //       std::vector<double>(shuffled_field.begin() + (i * prop_count),
+  //                           shuffled_field.begin() + ((i + 1) *
+  //                           prop_count)));
+  // }
 
-  this->setDumpedField(input);
-  PhreeqcRM::RunCells();
-  this->getDumpedField(input);
+  // this->setDumpedField(input);
+  // PhreeqcRM::RunCells();
+  // this->getDumpedField(input);
 
-  shuffled_field.clear();
-  for (std::size_t i = 0; i < n_cells; i++) {
-    shuffled_field.insert(shuffled_field.end(), input[i].begin(),
-                          input[i].end());
-  }
+  // shuffled_field.clear();
+  // for (std::size_t i = 0; i < n_cells; i++) {
+  //   shuffled_field.insert(shuffled_field.end(), input[i].begin(),
+  //                         input[i].end());
+  // }
 
-  std::vector<double> out_vec{shuffled_field};
-  unshuffleField(shuffled_field, n_cells, prop_count, 1, out_vec);
-  chem_field = out_vec;
+  // std::vector<double> out_vec{shuffled_field};
+  // unshuffleField(shuffled_field, n_cells, prop_count, 1, out_vec);
+  // chem_field = out_vec;
 }
 
-void poet::ChemistryModule::MasterRunParallel() {
+void poet::ChemistryModule::MasterRunParallel(double dt) {
   /* declare most of the needed variables here */
   double seq_a, seq_b, seq_c, seq_d;
   double worker_chemistry_a, worker_chemistry_b;
@@ -353,21 +368,30 @@ void poet::ChemistryModule::MasterRunParallel() {
   int pkg_to_send, pkg_to_recv;
   int free_workers;
   int i_pkgs;
+  int ftype;
 
-  int ftype = CHEM_WORK_LOOP;
+  const std::vector<uint32_t> wp_sizes_vector =
+      CalculateWPSizesVector(this->n_cells, this->wp_size);
+
+  if (this->ai_surrogate_enabled) {
+    ftype = CHEM_AI_BCAST_VALIDITY;
+    PropagateFunctionType(ftype);
+    this->ai_surrogate_validity_vector = shuffleVector(this->ai_surrogate_validity_vector,
+                                                       this->n_cells, 
+                                                       wp_sizes_vector.size());
+    ChemBCast(&this->ai_surrogate_validity_vector.front(), this->n_cells, MPI_INT);
+  }  
+
+  ftype = CHEM_WORK_LOOP;
   PropagateFunctionType(ftype);
 
   MPI_Barrier(this->group_comm);
 
-  double dt = this->PhreeqcRM::GetTimeStep();
   static uint32_t iteration = 0;
 
   /* start time measurement of sequential part */
   seq_a = MPI_Wtime();
 
-  const std::vector<uint32_t> wp_sizes_vector =
-      CalculateWPSizesVector(this->n_cells, this->wp_size);
-
   /* shuffle grid */
   // grid.shuffleAndExport(mpi_buffer);
   std::vector<double> mpi_buffer =
@@ -465,3 +489,13 @@ poet::ChemistryModule::CalculateWPSizesVector(uint32_t n_cells,
 
   return wp_sizes_vector;
 }
+
+void poet::ChemistryModule::masterSetField(Field field) {
+  this->chem_field = field;
+  this->prop_count = field.GetProps().size();
+
+  int ftype = CHEM_FIELD_INIT;
+  PropagateFunctionType(ftype);
+
+  ChemBCast(&this->prop_count, 1, MPI_UINT32_T);
+}

src/Chemistry/SurrogateModels/DHT.c

@@ -1,4 +1,3 @@
-/// Time-stamp: "Last modified 2023-08-10 11:50:46 mluebke"
 /*
 ** Copyright (C) 2017-2021 Max Luebke (University of Potsdam)
 **
@@ -16,8 +15,9 @@
 ** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */
 
+#include "DHT.h"
+
 #include <mpi.h>
-#include <poet/DHT.h>
 
 #include <inttypes.h>
 #include <math.h>
@@ -105,6 +105,8 @@ DHT *DHT_create(MPI_Comm comm, uint64_t size, unsigned int data_size,
   object->index_count = 9 - (index_bytes / 8);
   object->index = (uint64_t *)malloc((object->index_count) * sizeof(uint64_t));
   object->mem_alloc = mem_alloc;
+  object->sum_idx = 0;
+  object->cnt_idx = 0;
 
   // if set, initialize dht_stats
 #ifdef DHT_STATISTICS
@@ -189,6 +191,9 @@ int DHT_write_accumulate(DHT *table, const void *send_key, int data_size,
     }
   }
 
+  table->cnt_idx += 1;
+  table->sum_idx += (i + 1);
+
   if (result == DHT_WRITE_SUCCESS_WITH_EVICTION) {
     memset((char *)table->send_entry + 1 + table->key_size, '\0',
            table->data_size);
@@ -279,6 +284,9 @@ int DHT_write(DHT *table, void *send_key, void *send_data, uint32_t *proc,
     }
   }
 
+  table->cnt_idx += 1;
+  table->sum_idx += (i + 1);
+
   // put data to DHT (with last selected index by value i)
   if (MPI_Put(table->send_entry, 1 + table->data_size + table->key_size,
               MPI_BYTE, dest_rank, table->index[i],
@@ -550,6 +558,41 @@ int DHT_free(DHT *table, int *eviction_counter, int *readerror_counter) {
   return DHT_SUCCESS;
 }
 
+float DHT_get_used_idx_factor(DHT *table, int with_reset) {
+  int rank;
+  MPI_Comm_rank(table->communicator, &rank);
+
+  float my_avg_idx = (float)table->sum_idx / (float)table->cnt_idx;
+
+  float max_mean_index;
+
+  MPI_Reduce(&my_avg_idx, &max_mean_index, 1, MPI_FLOAT, MPI_MAX, 0,
+             table->communicator);
+
+  MPI_Bcast(&max_mean_index, 1, MPI_FLOAT, 0, table->communicator);
+
+  if (!!with_reset) {
+    table->sum_idx = 0;
+    table->cnt_idx = 0;
+  }
+
+  return max_mean_index;
+}
+
+int DHT_flush(DHT *table) {
+  // make sure all processes are synchronized
+  MPI_Barrier(table->communicator);
+
+  // wipe local memory with zeros
+  memset(table->mem_alloc, '\0',
+         table->table_size * (1 + table->data_size + table->key_size));
+
+  table->sum_idx = 0;
+  table->cnt_idx = 0;
+
+  return DHT_SUCCESS;
+}
+
 int DHT_print_statistics(DHT *table) {
 #ifdef DHT_STATISTICS
   int *written_buckets;

src/Chemistry/SurrogateModels/DHT.h

@@ -117,6 +117,9 @@ typedef struct {
   unsigned int index_count;
 
   int (*accumulate_callback)(int, void *, int, void *);
+
+  size_t sum_idx;
+  size_t cnt_idx;
 #ifdef DHT_STATISTICS
   /** Detailed statistics of the usage of the DHT. */
   DHT_stats *stats;
@@ -125,10 +128,11 @@ typedef struct {
 
 extern void DHT_set_accumulate_callback(DHT *table,
                                         int (*callback_func)(int, void *, int,
-                                                              void *));
+                                                             void *));
 
 extern int DHT_write_accumulate(DHT *table, const void *key, int send_size,
-                                void *data, uint32_t *proc, uint32_t *index, int *callback_ret);
+                                void *data, uint32_t *proc, uint32_t *index,
+                                int *callback_ret);
 
 /**
  * @brief Create a DHT.
@@ -284,44 +288,8 @@ extern int DHT_free(DHT *table, int *eviction_counter, int *readerror_counter);
  */
 extern int DHT_print_statistics(DHT *table);
 
-/**
- * @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 float DHT_get_used_idx_factor(DHT *table, int with_reset);
 
-/**
- * @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);
+extern int DHT_flush(DHT *table);
 
 #endif /* DHT_H */

src/Chemistry/SurrogateModels/DHT_Wrapper.cpp

@@ -1,4 +1,3 @@
-//  Time-stamp: "Last modified 2023-09-08 22:09:03 mluebke"
 
 /*
 ** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
@@ -20,22 +19,18 @@
 ** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */
 
-#include "poet/DHT_Wrapper.hpp"
-#include "poet/HashFunctions.hpp"
-#include "poet/Interpolation.hpp"
-#include "poet/LookupKey.hpp"
-#include "poet/Rounding.hpp"
-#include "poet/enums.hpp"
+#include "DHT_Wrapper.hpp"
 
-#include "poet/RInsidePOET.hpp"
+#include "Init/InitialList.hpp"
+#include "Rounding.hpp"
 
+#include <Rcpp/proxy/ProtectedProxy.h>
 #include <algorithm>
 #include <cassert>
 #include <cmath>
 #include <cstddef>
 #include <cstdint>
 #include <cstring>
-#include <iostream>
 #include <stdexcept>
 #include <vector>
 
@@ -47,12 +42,13 @@ DHT_Wrapper::DHT_Wrapper(MPI_Comm dht_comm, std::uint64_t dht_size,
                          const NamedVector<std::uint32_t> &key_species,
                          const std::vector<std::int32_t> &key_indices,
                          const std::vector<std::string> &_output_names,
-                         const ChemistryParams::Chem_Hook_Functions &_hooks,
-                         uint32_t data_count, bool _with_interp)
+                         const InitialList::ChemistryHookFunctions &_hooks,
+                         uint32_t data_count, bool _with_interp,
+                         bool _has_het_ids)
     : key_count(key_indices.size()), data_count(data_count),
       input_key_elements(key_indices), communicator(dht_comm),
       key_species(key_species), output_names(_output_names), hooks(_hooks),
-      with_interp(_with_interp) {
+      with_interp(_with_interp), has_het_ids(_has_het_ids) {
   // initialize DHT object
   // key size = count of key elements + timestep
   uint32_t key_size = (key_count + 1) * sizeof(Lookup_Keyelement);
@@ -133,43 +129,47 @@ void DHT_Wrapper::fillDHT(const WorkPackage &work_package) {
   dht_results.filledDHT = std::vector<bool>(length, false);
   for (int i = 0; i < length; i++) {
     // If true grid cell was simulated, needs to be inserted into dht
-    if (work_package.mapping[i] == CHEM_PQC) {
-
-      // check if calcite or dolomite is absent and present, resp.n and vice
-      // versa in input/output. If this is the case -> Do not write to DHT!
-      // HACK: hardcoded, should be fixed!
-      if (hooks.dht_fill.isValid()) {
-        NamedVector<double> old_values(output_names, work_package.input[i]);
-        NamedVector<double> new_values(output_names, work_package.output[i]);
-
-        if (hooks.dht_fill(old_values, new_values)) {
-          continue;
-        }
-      }
-
-      uint32_t proc, index;
-      auto &key = dht_results.keys[i];
-      const auto data =
-          (with_interp
-               ? outputToInputAndRates(work_package.input[i],
-                                       work_package.output[i])
-               : work_package.output[i]);
-      // void *data = (void *)&(work_package[i * this->data_count]);
-      // fuzz data (round, logarithm etc.)
-
-      // insert simulated data with fuzzed key into DHT
-      int res = DHT_write(this->dht_object, key.data(),
-                          const_cast<double *>(data.data()), &proc, &index);
-
-      dht_results.locations[i] = {proc, index};
-
-      // if data was successfully written ...
-      if ((res != DHT_SUCCESS) && (res == DHT_WRITE_SUCCESS_WITH_EVICTION)) {
-        dht_evictions++;
-      }
-
-      dht_results.filledDHT[i] = true;
+    if (work_package.mapping[i] != CHEM_PQC) {
+      continue;
     }
+
+    if (work_package.input[i][0] != 2) {
+      continue;
+    }
+
+    // check if calcite or dolomite is absent and present, resp.n and vice
+    // versa in input/output. If this is the case -> Do not write to DHT!
+    // HACK: hardcoded, should be fixed!
+    if (hooks.dht_fill.isValid()) {
+      NamedVector<double> old_values(output_names, work_package.input[i]);
+      NamedVector<double> new_values(output_names, work_package.output[i]);
+
+      if (Rcpp::as<bool>(hooks.dht_fill(old_values, new_values))) {
+        continue;
+      }
+    }
+
+    uint32_t proc, index;
+    auto &key = dht_results.keys[i];
+    const auto data =
+        (with_interp ? outputToInputAndRates(work_package.input[i],
+                                             work_package.output[i])
+                     : work_package.output[i]);
+    // void *data = (void *)&(work_package[i * this->data_count]);
+    // fuzz data (round, logarithm etc.)
+
+    // insert simulated data with fuzzed key into DHT
+    int res = DHT_write(this->dht_object, key.data(),
+                        const_cast<double *>(data.data()), &proc, &index);
+
+    dht_results.locations[i] = {proc, index};
+
+    // if data was successfully written ...
+    if ((res != DHT_SUCCESS) && (res == DHT_WRITE_SUCCESS_WITH_EVICTION)) {
+      dht_evictions++;
+    }
+
+    dht_results.filledDHT[i] = true;
   }
 }
 
@@ -273,9 +273,10 @@ LookupKey DHT_Wrapper::fuzzForDHT_R(const std::vector<double> &cell,
 
   NamedVector<double> input_nv(this->output_names, cell);
 
-  const std::vector<double> eval_vec = hooks.dht_fuzz(input_nv);
+  const std::vector<double> eval_vec =
+      Rcpp::as<std::vector<double>>(hooks.dht_fuzz(input_nv));
   assert(eval_vec.size() == this->key_count);
-  LookupKey vecFuzz(this->key_count + 1, {.0});
+  LookupKey vecFuzz(this->key_count + 1 + has_het_ids, {.0});
 
   DHT_Rounder rounder;
 
@@ -295,6 +296,9 @@ LookupKey DHT_Wrapper::fuzzForDHT_R(const std::vector<double> &cell,
   }
   // add timestep to the end of the key as double value
   vecFuzz[this->key_count].fp_element = dt;
+  if (has_het_ids) {
+    vecFuzz[this->key_count + 1].fp_element = cell[0];
+  }
 
   return vecFuzz;
 }
@@ -302,7 +306,7 @@ LookupKey DHT_Wrapper::fuzzForDHT_R(const std::vector<double> &cell,
 LookupKey DHT_Wrapper::fuzzForDHT(const std::vector<double> &cell, double dt) {
   const auto c_zero_val = std::pow(10, AQUEOUS_EXP);
 
-  LookupKey vecFuzz(this->key_count + 1, {.0});
+  LookupKey vecFuzz(this->key_count + 1 + has_het_ids, {.0});
   DHT_Rounder rounder;
 
   int totals_i = 0;
@@ -328,6 +332,9 @@ LookupKey DHT_Wrapper::fuzzForDHT(const std::vector<double> &cell, double dt) {
   }
   // add timestep to the end of the key as double value
   vecFuzz[this->key_count].fp_element = dt;
+  if (has_het_ids) {
+    vecFuzz[this->key_count + 1].fp_element = cell[0];
+  }
 
   return vecFuzz;
 }

src/Chemistry/SurrogateModels/DHT_Wrapper.hpp

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

src/Chemistry/SurrogateModels/HashFunctions.cpp

@@ -1,4 +1,3 @@
-//  Time-stamp: "Last modified 2023-04-24 23:20:55 mluebke"
 /*
 **-----------------------------------------------------------------------------
 ** MurmurHash2 was written by Austin Appleby, and is placed in the public
@@ -24,7 +23,7 @@
 ** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */
 
-#include "poet/HashFunctions.hpp"
+#include "HashFunctions.hpp"
 
 #if defined(_MSC_VER)
 

src/Chemistry/SurrogateModels/HashFunctions.hpp

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

src/Chemistry/SurrogateModels/Interpolation.hpp

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