Merge branch 'mdl/defaultstorelease' into 'main'

DEFAULT_BUILD_TYPE set to Release

See merge request naaice/poet!60

.devcontainer/Dockerfile

@@ -1,20 +1,108 @@
-FROM mcr.microsoft.com/vscode/devcontainers/base:debian
+FROM gcc:11.2.0 AS builder
 
-RUN sudo apt-get update && export DEBIAN_FRONTEND=noninteractive \
-    && sudo apt-get install -y \
-        cmake \
-        git \
-        libeigen3-dev \
-        libopenmpi-dev \
-        r-base-dev 
+ENV DEBIAN_FRONTEND=noninteractive
 
-RUN git clone https://github.com/doctest/doctest.git /doctest \
-    && cd /doctest \
-    && mkdir build \
-    && cd build \
-    && cmake .. \
-    && make install \
-    && cd / \
-    && rm -rf /doctest
+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/*
 
-RUN  /usr/bin/R -q -e "install.packages(c('Rcpp', 'RInside'))"
+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

@@ -21,7 +21,8 @@
 	},
 	// 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}/.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"

.gitlab-ci.yml

@@ -20,7 +20,6 @@ image: git.gfz-potsdam.de:5000/naaice/poet:ci
 
 stages:          # List of stages for jobs, and their order of execution
   - release
-  - build
   - test
 
 variables:
@@ -28,26 +27,11 @@ 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
+test:       # This job runs in the build stage, which runs first.
+  stage: test
   script:
     - mkdir -p build && cd build
-    - cmake -DPOET_ENABLE_TESTING=ON -DPOET_PREPROCESS_BENCHS=OFF ..
-    - make -j$(nproc)
-  cache: 
-    key: build-$CI_COMMIT_BRANCH
-    paths:
-      - build  
-  artifacts:
-    paths:
-      - build
-
-test-poet:
-  stage: test
-  dependencies:
-    - build-poet
-  script:
-    - cd build
+    - cmake -DPOET_ENABLE_TESTING=ON -DPOET_PREPROCESS_BENCHS=OFF -DCMAKE_BUILD_TYPE=Release ..
     - make -j$(nproc) check
 
 pages:
@@ -65,6 +49,22 @@ pages:
   rules:
     - if: $CI_COMMIT_REF_NAME == $CI_DEFAULT_BRANCH || $CI_COMMIT_TAG
 
+push:
+  stage: release
+  variables:
+    GITHUB_REPOSITORY: 'git@github.com:POET-Simulator/POET.git'
+  before_script:
+    # I know that there is this file env variable in gitlab, but somehow it does not work for me (still complaining about white spaces ...)
+    # Therefore, the ssh key is stored as a base64 encoded string
+    - mkdir -p ~/.ssh && echo $GITHUB_SSH_PRIVATE_KEY | base64 -d > ~/.ssh/id_ed25519 && chmod 0600 ~/.ssh/id_ed25519
+    - ssh-keyscan github.com >> ~/.ssh/known_hosts
+    - echo $MIRROR_SCRIPT | base64 -d > mirror.sh && chmod +x mirror.sh
+  script:
+    - if [[-d poet.git ]]; then rm -rf poet.git; fi
+    - git clone --mirror "https://git.gfz-potsdam.de/naaice/poet.git" "poet.git" && cd poet.git
+    - git push --mirror $GITHUB_REPOSITORY
+  allow_failure: true
+
 #archive-sources:       # This job runs in the build stage, which runs first.
 #  image: python:3
 #  stage: release

.gitmodules

@@ -2,6 +2,6 @@
 	path = ext/tug
 	url = ../tug.git
 
-[submodule "ext/iphreeqc"]
-	path = ext/iphreeqc
-	url = ../iphreeqc.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,4 +1,4 @@
-cmake_minimum_required(VERSION 3.14)
+cmake_minimum_required(VERSION 3.20)
 
 project(POET
   LANGUAGES CXX C
@@ -8,6 +8,16 @@ project(POET
 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")
@@ -31,7 +41,7 @@ endif()
 set(TUG_ENABLE_TESTING OFF CACHE BOOL "" FORCE)
 
 add_subdirectory(ext/tug EXCLUDE_FROM_ALL)
-add_subdirectory(ext/iphreeqc EXCLUDE_FROM_ALL)
+add_subdirectory(ext/litephreeqc EXCLUDE_FROM_ALL)
 
 option(POET_ENABLE_TESTING "Build test suite for POET" OFF)
 
@@ -41,4 +51,4 @@ endif()
 
 option(BUILD_DOC "Build documentation with doxygen" OFF)
 
-add_subdirectory(docs)
+add_subdirectory(docs)

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,80 +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)
-- **IPhreeqc** with patches from GFZ -
-  https://github.com/usgs-coupled/iphreeqc -
-  https://git.gfz-potsdam.de/naaice/iphreeqc
-- **tug** - https://git.gfz-potsdam.de/naaice/tug
+- **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+
+- CMake 3.20+
 - Eigen3 3.4+ (required by `tug`)
-- *optional*: `doxygen` with `dot` bindings for documentiation
+- *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_PREPROCESS_BENCHS**=*boolean* - enables the preprocessing of predefined
-  models/benchmarks. Defaults to *ON*.
+- **POET_DHT_Debug**=_boolean_ - toggles the output of detailed
+  statistics about DHT usage. Defaults to _OFF_.
+- **POET_ENABLE_TESTING**=_boolean_ - enables small set of unit tests
+  (more to come). Defaults to _OFF_.
+- **POET_PHT_ADDITIONAL_INFO**=_boolean_ - enabling the count of
+  accesses to one PHT bucket. Use with caution, as things will get
+  slowed down significantly. Defaults to _OFF_.
+- **POET_PREPROCESS_BENCHS**=*boolean* - enables the preprocessing of
+  predefined models/benchmarks. Defaults to *ON*.
   
 ### Example: Build from scratch
 
@@ -87,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
@@ -95,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
 ```
@@ -115,53 +143,59 @@ poet
 └── share
     └── poet
         ├── barite
-        │   ├── barite_200.rds
+        │   ├── barite_200.qs2
         │   ├── barite_200_rt.R
-        │   ├── barite_het.rds
+        │   ├── barite_het.qs2
         │   └── barite_het_rt.R
         ├── dolo
-        │   ├── dolo_inner_large.rds
+        │   ├── dolo_inner_large.qs2
         │   ├── dolo_inner_large_rt.R
-        │   ├── dolo_interp.rds
+        │   ├── dolo_interp.qs2
         │   └── dolo_interp_rt.R
         └── surfex
-            ├── PoetEGU_surfex_500.rds
+            ├── PoetEGU_surfex_500.qs2
             └── PoetEGU_surfex_500_rt.R
 ```
 
 With the installation of POET, two executables are provided: 
   - `poet` - the main executable to run simulations
-  - `poet_init` - a preprocessor to generate input files for POET from R scripts
+  - `poet_init` - a preprocessor to generate input files for POET from
+    R scripts
 
-Preprocessed benchmarks can be found in the `share/poet` directory with an
-according *runtime* setup. More on those files and how to create them later. 
+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] <RUNFILE> <SIMFILE> <OUTPUT_DIRECTORY>`
-where:
+Run POET by `mpirun ./poet [OPTIONS] <RUNFILE> <SIMFILE>
+<OUTPUT_DIRECTORY>` where:
 
 - **OPTIONS** - POET options (explained below)
-- **RUNFILE** - Runtime parameters described as R script 
+- **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
+- **OUTPUT_DIRECTORY** - path, where all output of POET should be
+  stored
 
-### POET 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_)                                                                                  |
-| **-P, --progress**          |              | show progress bar                                                                                                        |
-| **--dht**                   |              | enabling DHT usage (defaults to _OFF_)                                                                                   |
-| **--dht-strategy=**         | _0-1_        | change DHT strategy. **NOT IMPLEMENTED YET** (Defaults to _0_)                                                           |
-| **--dht-size=**             | _1-n_        | size of DHT per process involved in megabyte (defaults to _1000 MByte_)                                                  |
-| **--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`
 
@@ -175,39 +209,91 @@ Following values can be set:
 
 ### Example: Running from scratch
 
-We will continue the above example and start a simulation with *barite_het*,
-which simulation files can be found in
-`<INSTALL_DIR>/share/poet/barite/barite_het*`. As transport a heterogeneous
-diffusion is used. It's a small 2D grid, 2x5 grid, simulating 50 time steps with
-a time step size of 100 seconds. To start the simulation with 4 processes `cd`
-into your previously installed POET-dir `<POET_INSTALL_DIR>/bin` and run:
+We will continue the above example and start a simulation with
+*barite_het*, which simulation files can be found in
+`<INSTALL_DIR>/share/poet/barite/barite_het*`. As transport a
+heterogeneous diffusion is used. It's a small 2D grid, 2x5 grid,
+simulating 50 time steps with a time step size of 100 seconds. To
+start the simulation with 4 processes `cd` into your previously
+installed POET-dir `<POET_INSTALL_DIR>/bin` and run:
 
 ```sh
 cp ../share/poet/barite/barite_het* .
-mpirun -n 4 ./poet barite_het_rt.R barite_het.rds output
+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 barite_het_rt.R barite_het.rds 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-potsdam.de/naaice/poet/-/wikis/Initialization). We try to
-keep the document up-to-date. However, if you encounter missing information or
-need help, please get in touch with us via the issue tracker or E-Mail.
+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:
 
@@ -217,11 +303,12 @@ need help, please get in touch with us via the issue tracker or E-Mail.
 
 where: 
 
-- **output** - name of the output file (defaults to the input file name
-  with the extension `.rds`)
-- **setwd** - set the working directory to the directory of the input file (e.g.
-  to allow relative paths in the input script). However, the output file
-  will be stored in the directory from which `poet_init` was called.
+- **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()
 
@@ -231,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/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

@@ -1,21 +1,67 @@
-pqc_to_grid <- function(pqc_in, grid) {
+### 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_in <- as.matrix(pqc_in)
+    pqc_mat <- as.matrix(pqc_mat)
 
     # Flatten the matrix into a vector
-    id_vector <- as.numeric(t(grid))
+    id_vector <- as.integer(t(grid))
 
     # Find the matching rows in the matrix
-    row_indices <- match(id_vector, pqc_in[, "ID"])
+    row_indices <- match(id_vector, pqc_mat[, "ID"])
 
-    # Extract the matching rows from pqc_in to size of grid matrix
-    result_mat <- pqc_in[row_indices, ]
+    # 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)]
+    # res_df <- res_df[, colSums(is.na(res_df)) != nrow(res_df)]
 
     # Remove row names
     rownames(res_df) <- NULL
@@ -23,6 +69,12 @@ pqc_to_grid <- function(pqc_in, grid) {
     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]
@@ -34,6 +86,10 @@ resolve_pqc_bound <- function(pqc_mat, transport_spec, id) {
     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)
@@ -53,4 +109,4 @@ add_missing_transport_species <- function(init_grid, 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
@@ -35,14 +33,18 @@ master_init <- function(setup, out_dir, init_field) {
     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_0.rds")
+        init_field_out <- paste0(out_dir, "/iter_", sprintf(fmt = fmt, 0), ".", setup$out_ext)
         init_field <- data.frame(init_field, check.names = FALSE)
-        saveRDS(init_field, file = init_field_out)
+        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)
@@ -61,7 +63,7 @@ master_iteration_end <- function(setup, state_T, state_C) {
     iter <- setup$iter
     # print(iter)
     ## max digits for iterations
-    dgts <- as.integer(ceiling(log10(setup$maxiter)))
+    dgts <- as.integer(ceiling(log10(setup$maxiter + 1)))
     ## string format to use in sprintf
     fmt <- paste0("%0", dgts, "d")
 
@@ -69,23 +71,30 @@ master_iteration_end <- function(setup, state_T, state_C) {
     ## comprised in setup$out_save
     if (setup$store_result) {
         if (iter %in% setup$out_save) {
-            nameout <- paste0(setup$out_dir, "/iter_", sprintf(fmt = fmt, iter), ".rds")
-
+            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)
 
-            saveRDS(list(
+            ai_surrogate_info <- list(
+                prediction_time = if (exists("ai_prediction_time")) as.integer(ai_prediction_time) else NULL,
+                training_time = if (exists("ai_training_time")) as.integer(ai_training_time) else NULL,
+                valid_predictions = if (exists("validity_vector")) validity_vector else NULL
+            )
+
+            SaveRObj(x = list(
                 T = state_T,
                 C = state_C,
-                simtime = as.integer(setup$simulation_time)
-            ), file = nameout)
+                simtime = as.integer(setup$simulation_time),
+                totaltime = as.integer(totaltime),
+                ai_surrogate_info = ai_surrogate_info
+            ), path = nameout)
             msgm("results stored in <", nameout, ">")
         }
     }
     ## Add last time step to simulation time
     setup$simulation_time <- setup$simulation_time + setup$timesteps[iter]
 
-    msgm("done iteration", iter, "/", length(setup$timesteps))
+    ## msgm("done iteration", iter, "/", length(setup$timesteps))
     setup$iter <- setup$iter + 1
     return(setup)
 }
@@ -100,70 +109,103 @@ msgm <- function(...) {
 }
 
 
-## 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"))
-}
-
 GetWorkPackageSizesVector <- function(n_packages, package_size, 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"))
+## }

bench/CMakeLists.txt

@@ -1,4 +1,3 @@
-
 function(ADD_BENCH_TARGET TARGET POET_BENCH_LIST RT_FILES OUT_PATH)
     set(bench_install_dir share/poet/${OUT_PATH})
 
@@ -7,18 +6,19 @@ function(ADD_BENCH_TARGET TARGET POET_BENCH_LIST RT_FILES OUT_PATH)
 
     foreach(BENCH_FILE ${${POET_BENCH_LIST}})
         get_filename_component(BENCH_NAME ${BENCH_FILE} NAME_WE)
-        set(OUT_FILE ${CMAKE_CURRENT_BINARY_DIR}/${BENCH_NAME}.rds)
+        set(OUT_FILE ${CMAKE_CURRENT_BINARY_DIR}/${BENCH_NAME})
+        set(OUT_FILE_EXT ${OUT_FILE}.qs2)
 
         add_custom_command(
-            OUTPUT ${OUT_FILE}
-            COMMAND $<TARGET_FILE:poet_init> -o ${OUT_FILE} -s ${CMAKE_CURRENT_SOURCE_DIR}/${BENCH_FILE}
+            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})
+        list(APPEND OUT_FILES_LIST ${OUT_FILE_EXT})
 
     endforeach(BENCH_FILE ${${POET_BENCH_LIST}})
         
@@ -41,4 +41,3 @@ add_custom_target(${BENCHTARGET} ALL)
 add_subdirectory(barite)
 add_subdirectory(dolo)
 add_subdirectory(surfex)
-

bench/barite/README.org

@@ -22,6 +22,15 @@ mpirun -np 4 ./poet --interp  barite_interp_eval.R barite_results
   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_200.R

@@ -35,19 +35,20 @@ diffusion_setup <- list(
 )
 
 dht_species <- c(
-  "H"         = 7,
-  "O"         = 7,
-  "Charge"    = 4,
-  "Ba"        = 7,
-  "Cl"        = 7,
-  "S(6)"      = 7,
-  "Sr"        = 7,
-  "Barite"    = 4,
-  "Celestite" = 4
+  "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
+  dht_species = dht_species,
+  ai_surrogate_input_script = "./barite_200ai_surrogate_input_script.R"
 )
 
 # Define a setup list for simulation configuration

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/dolo/dolo_interp.R

@@ -41,36 +41,36 @@ diffusion_setup <- list(
 )
 
 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)
-    }
+    # 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)])
-}
+# 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(4)" = 6,
+        "C" = 6,
         "Ca" = 6,
         "Cl" = 3,
         "Mg" = 5,
@@ -95,7 +95,7 @@ check_neg_cal_dol <- function(result) {
 # significant digits)
 
 pht_species <- c(
-    "C(4)" = 3,
+    "C" = 3,
     "Ca" = 3,
     "Mg" = 2,
     "Calcite" = 2,
@@ -107,7 +107,7 @@ chemistry_setup <- list(
         "H" = 3,
         "O" = 3,
         "Charge" = 3,
-        "C(4)" = 6,
+        "C" = 6,
         "Ca" = 6,
         "Cl" = 3,
         "Mg" = 5,
@@ -117,7 +117,7 @@ chemistry_setup <- list(
     pht_species = pht_species,
     hooks = list(
         dht_fill = check_sign_cal_dol_dht,
-        dht_fuzz = fuzz_input_dht_keys,
+        # dht_fuzz = fuzz_input_dht_keys,
         interp_pre = check_sign_cal_dol_interp,
         interp_post = check_neg_cal_dol
     )

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

ext/iphreeqc

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

ext/litephreeqc

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

ext/tug

@@ -1 +1 @@
-Subproject commit 449647010ab9cdf9e405139f360424a2b21ab3ab
+Subproject commit 9c4aeee410c71d064f7567143d4f8d6451ade75a

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

@@ -1,17 +1,13 @@
-#ifndef RPOET_H_
-#define RPOET_H_
+#pragma once
 
 #include <RInside.h>
 #include <Rcpp.h>
 #include <Rinternals.h>
-#include <cstddef>
 #include <exception>
-#include <optional>
-#include <stdexcept>
+#include <memory>
 #include <string>
-#include <utility>
-#include <vector>
 
+namespace poet {
 class RInsidePOET : public RInside {
 public:
   static RInsidePOET &getInstance() {
@@ -33,44 +29,64 @@ private:
   RInsidePOET() : RInside(){};
 };
 
-template <typename T> class RHookFunction {
+/**
+ * @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:
-  RHookFunction() {}
-  RHookFunction(RInside &R, const std::string &f_name) {
+  DEFunc() {}
+  DEFunc(const std::string &f_name) {
     try {
-      this->func = Rcpp::Function(Rcpp::as<SEXP>(R.parseEval(f_name.c_str())));
+      this->func = std::make_shared<Rcpp::Function>(f_name);
     } catch (const std::exception &e) {
     }
   }
 
-  RHookFunction(SEXP f) {
+  DEFunc(SEXP f) {
     try {
-      this->func = Rcpp::Function(f);
+      this->func = std::make_shared<Rcpp::Function>(f);
     } 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...)));
+  template <typename... Args> SEXP operator()(Args... args) const {
+    if (func) {
+      return (*this->func)(args...);
     } else {
       throw std::exception();
     }
   }
 
-  RHookFunction &operator=(const RHookFunction &rhs) {
+  DEFunc &operator=(const DEFunc &rhs) {
     this->func = rhs.func;
     return *this;
   }
 
-  RHookFunction(const RHookFunction &rhs) { this->func = rhs.func; }
+  DEFunc(const DEFunc &rhs) { this->func = rhs.func; }
 
-  bool isValid() const { return this->func.has_value(); }
+  bool isValid() const { return static_cast<bool>(func); }
 
-  SEXP asSEXP() const { return Rcpp::as<SEXP>(this->func.value()); }
+  SEXP asSEXP() const {
+    if (!func) {
+      return R_NilValue;
+    }
+    return Rcpp::as<SEXP>(*this->func.get());
+  }
 
 private:
-  std::optional<Rcpp::Function> func;
+  std::shared_ptr<Rcpp::Function> func;
 };
 
-#endif // RPOET_H_
+} // namespace poet

src/Base/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));
-   }
-}
-
-

src/CMakeLists.txt

@@ -1,7 +1,4 @@
-add_library(POETLib)
-
-target_sources(POETLib 
-  PRIVATE
+add_library(POETLib
   Init/InitialList.cpp 
   Init/GridInit.cpp 
   Init/DiffusionInit.cpp
@@ -18,6 +15,15 @@ target_sources(POETLib
   Chemistry/SurrogateModels/ProximityHashTable.cpp
 )
 
+set(POET_TUG_APPROACH "Implicit" CACHE STRING "tug numerical approach to use")
+set_property(CACHE POET_TUG_APPROACH PROPERTY STRINGS "Implicit" "Explicit")
+
+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()
+
 target_include_directories(POETLib PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}")
 target_link_libraries(
   POETLib 
@@ -27,6 +33,16 @@ target_link_libraries(
   PUBLIC MPI::MPI_C
 )
 
+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
@@ -70,16 +86,16 @@ target_compile_definitions(POETLib PUBLIC STRICT_R_HEADERS OMPI_SKIP_MPICXX)
 
 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)
+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)
+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

@@ -6,7 +6,7 @@
 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 };
+enum CHEMISTRY_OUT_SOURCE { CHEM_PQC, CHEM_DHT, CHEM_INTERP, CHEM_AISURR };
 
 struct WorkPackage {
   std::size_t size;
@@ -14,10 +14,7 @@ struct WorkPackage {
   std::vector<std::vector<double>> output;
   std::vector<std::uint8_t> mapping;
 
-  WorkPackage(std::size_t _size) : size(_size) {
-    input.resize(size);
-    output.resize(size);
-    mapping.resize(size, CHEM_PQC);
-  }
+  WorkPackage(std::size_t _size)
+      : size(_size), input(size), output(size), mapping(size, CHEM_PQC) {}
 };
 } // namespace poet

src/Chemistry/ChemistryModule.cpp

@@ -1,11 +1,14 @@
 #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>
@@ -63,7 +66,8 @@ inverseDistanceWeighting(const std::vector<std::int32_t> &to_calc,
       distance += std::pow(
           rescaled[key_comp_i][point_i] - rescaled[key_comp_i][data_set_n], 2);
     }
-    weights[point_i] = 1 / std::sqrt(distance);
+
+    weights[point_i] = distance != 0 ? 1 / std::sqrt(distance) : 0;
     assert(!std::isnan(weights[point_i]));
     inv_sum += weights[point_i];
   }
@@ -94,63 +98,9 @@ inverseDistanceWeighting(const std::vector<std::int32_t> &to_calc,
     key_delta /= inv_sum;
 
     results[output_comp_i] = from[output_comp_i] + key_delta;
+    assert(!std::isnan(results[output_comp_i]));
   }
 
-  // if (!has_h) {
-  //   double new_val = 0;
-  //   for (int j = 0; j < data_set_n; j++) {
-  //     new_val += weights[j] * output[j][0];
-  //   }
-  //   results[0] = new_val / inv_sum;
-  // }
-
-  // if (!has_h) {
-  //   double new_val = 0;
-  //   for (int j = 0; j < data_set_n; j++) {
-  //     new_val += weights[j] * output[j][1];
-  //   }
-  //   results[1] = new_val / inv_sum;
-  // }
-
-  // for (std::uint32_t i = 0; i < to_calc.size(); i++) {
-  //   const std::uint32_t interp_i = to_calc[i];
-
-  //   // rescale input between 0 and 1
-  //   for (int j = 0; j < input.size(); j++) {
-  //     buffer[j] = input[j].at(i);
-  //   }
-
-  //   buffer[buffer_size - 1] = from[interp_i];
-
-  //   const double min = *std::min_element(buffer, buffer + buffer_size);
-  //   const double max = *std::max_element(buffer, buffer + buffer_size);
-
-  //   for (int j = 0; j < input.size(); j++) {
-  //     buffer[j] = ((max - min) != 0 ? (buffer[j] - min) / (max - min) : 1);
-  //   }
-  //   from_rescaled =
-  //       ((max - min) != 0 ? (from[interp_i] - min) / (max - min) : 0);
-
-  //   double inv_sum = 0;
-
-  //   // calculate distances for each point
-  //   for (int i = 0; i < input.size(); i++) {
-  //     const double distance = std::pow(buffer[i] - from_rescaled, 2);
-
-  //     buffer[i] = distance > 0 ? (1 / std::sqrt(distance)) : 0;
-  //     inv_sum += buffer[i];
-  //   }
-  //   // calculate new values
-  //   double new_val = 0;
-  //   for (int i = 0; i < output.size(); i++) {
-  //     new_val += buffer[i] * output[i][interp_i];
-  //   }
-  //   results[interp_i] = new_val / inv_sum;
-  //   if (std::isnan(results[interp_i])) {
-  //     std::cout << "nan with new_val = " << output[0][i] << std::endl;
-  //   }
-  // }
-
   return results;
 }
 
@@ -167,16 +117,12 @@ poet::ChemistryModule::ChemistryModule(
   this->n_cells = chem_params.total_grid_cells;
 
   if (!is_master) {
-    for (const auto &[pqc_id, pqc_config] : chem_params.pqc_config) {
-      this->phreeqc_instances[pqc_id] =
-          std::make_unique<PhreeqcEngine>(pqc_config);
-    }
-    // for (std::size_t i = 0; i < chem_params.pqc_ids.size(); i++) {
-    //   this->phreeqc_instances[chem_params.pqc_ids[i]] =
-    //       std::make_unique<PhreeqcWrapper>(
-    //           chem_params.database, chem_params.pqc_scripts[i],
-    //           chem_params.pqc_sol_order, chem_params.field_header, wp_size_);
-    // }
+    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));
   }
 }
 
@@ -187,11 +133,10 @@ poet::ChemistryModule::~ChemistryModule() {
 }
 
 void poet::ChemistryModule::initializeDHT(
-    uint32_t size_mb, const NamedVector<std::uint32_t> &key_species) {
+    uint32_t size_mb, const NamedVector<std::uint32_t> &key_species,
+    bool has_het_ids) {
   constexpr uint32_t MB_FACTOR = 1E6;
 
-  this->dht_enabled = true;
-
   MPI_Comm dht_comm;
 
   if (this->is_master) {
@@ -221,7 +166,7 @@ void poet::ChemistryModule::initializeDHT(
 
     this->dht = new DHT_Wrapper(dht_comm, dht_size, map_copy, key_indices,
                                 this->prop_names, params.hooks,
-                                this->prop_count, interp_enabled);
+                                this->prop_count, interp_enabled, has_het_ids);
     this->dht->setBaseTotals(base_totals.at(0), base_totals.at(1));
   }
 }
@@ -312,9 +257,10 @@ void poet::ChemistryModule::initializeInterp(
       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);
+            static_cast<std::uint32_t>(map_copy[i]) -
+            (map_copy[i] > InterpolationModule::COARSE_DIFF
+                 ? InterpolationModule::COARSE_DIFF
+                 : 0);
         map_copy[i] = signif;
       }
     }
@@ -371,3 +317,8 @@ void poet::ChemistryModule::unshuffleField(const std::vector<double> &in_buffer,
     }
   }
 }
+
+void poet::ChemistryModule::set_ai_surrogate_validity_vector(
+    std::vector<int> r_vector) {
+  this->ai_surrogate_validity_vector = r_vector;
+}

src/Chemistry/ChemistryModule.hpp

@@ -8,10 +8,11 @@
 #include "ChemistryDefs.hpp"
 
 #include "Init/InitialList.hpp"
+#include "NameDouble.h"
 #include "SurrogateModels/DHT_Wrapper.hpp"
 #include "SurrogateModels/Interpolation.hpp"
 
-#include "PhreeqcEngine.hpp"
+#include "PhreeqcRunner.hpp"
 #include <array>
 #include <cstdint>
 #include <map>
@@ -75,14 +76,19 @@ public:
 
   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) {
@@ -92,9 +98,17 @@ public:
 
     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);
+      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) {
@@ -219,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;
@@ -228,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);
 
@@ -249,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 };
@@ -348,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 {
@@ -378,7 +403,7 @@ protected:
 
   const InitialList::ChemistryInit params;
 
-  std::map<int, std::unique_ptr<PhreeqcEngine>> phreeqc_instances;
+  std::unique_ptr<PhreeqcRunner> pqc_runner;
 };
 } // namespace poet
 

src/Chemistry/MasterFunctions.cpp

@@ -159,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,
@@ -247,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,
@@ -352,8 +368,21 @@ void poet::ChemistryModule::MasterRunParallel(double dt) {
   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);
@@ -363,9 +392,6 @@ void poet::ChemistryModule::MasterRunParallel(double dt) {
   /* 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 =

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)
 **
@@ -106,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
@@ -190,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);
@@ -280,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],
@@ -551,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-11-01 10:54:45 mluebke"
 
 /*
 ** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
@@ -25,6 +24,7 @@
 #include "Init/InitialList.hpp"
 #include "Rounding.hpp"
 
+#include <Rcpp/proxy/ProtectedProxy.h>
 #include <algorithm>
 #include <cassert>
 #include <cmath>
@@ -43,11 +43,12 @@ DHT_Wrapper::DHT_Wrapper(MPI_Comm dht_comm, std::uint64_t dht_size,
                          const std::vector<std::int32_t> &key_indices,
                          const std::vector<std::string> &_output_names,
                          const InitialList::ChemistryHookFunctions &_hooks,
-                         uint32_t data_count, bool _with_interp)
+                         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);
@@ -128,42 +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;
   }
 }
 
@@ -267,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;
 
@@ -289,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;
 }
@@ -296,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;
@@ -322,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
@@ -88,7 +87,7 @@ public:
               const std::vector<std::int32_t> &key_indices,
               const std::vector<std::string> &output_names,
               const InitialList::ChemistryHookFunctions &hooks,
-              uint32_t data_count, bool with_interp);
+              uint32_t data_count, bool with_interp, bool has_het_ids);
   /**
    * @brief Destroy the dht wrapper object
    *
@@ -265,6 +264,7 @@ private:
   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

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,4 +1,3 @@
-//  Time-stamp: "Last modified 2023-08-16 16:49:31 mluebke"
 
 #ifndef INTERPOLATION_H_
 #define INTERPOLATION_H_
@@ -167,6 +166,8 @@ public:
 
   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,6 +263,8 @@ private:
   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
 

src/Chemistry/SurrogateModels/InterpolationModule.cpp

@@ -1,4 +1,3 @@
-//  Time-stamp: "Last modified 2023-08-16 17:02:31 mluebke"
 #include "Init/InitialList.hpp"
 #include "Interpolation.hpp"
 
@@ -9,11 +8,13 @@
 #include "Rounding.hpp"
 
 #include <Rcpp.h>
+#include <Rcpp/proxy/ProtectedProxy.h>
 #include <Rinternals.h>
 
 #include <algorithm>
 #include <array>
 #include <cassert>
+#include <cmath>
 #include <cstddef>
 #include <cstdint>
 #include <functional>
@@ -75,6 +76,11 @@ void InterpolationModule::tryInterpolation(WorkPackage &work_package) {
   const auto dht_results = this->dht_instance.getDHTResults();
 
   for (int wp_i = 0; wp_i < work_package.size; wp_i++) {
+    if (work_package.input[wp_i][0] != 2) {
+      interp_result.status[wp_i] = INSUFFICIENT_DATA;
+      continue;
+    }
+
     if (work_package.mapping[wp_i] != CHEM_PQC) {
       interp_result.status[wp_i] = NOT_NEEDED;
       continue;
@@ -94,7 +100,8 @@ void InterpolationModule::tryInterpolation(WorkPackage &work_package) {
     if (hooks.interp_pre.isValid()) {
       NamedVector<double> nv_in(this->out_names, work_package.input[wp_i]);
 
-      auto rm_indices = hooks.interp_pre(nv_in, pht_result.in_values);
+      std::vector<int> rm_indices = Rcpp::as<std::vector<int>>(
+          hooks.interp_pre(nv_in, pht_result.in_values));
 
       pht_result.size -= rm_indices.size();
 
@@ -115,15 +122,30 @@ void InterpolationModule::tryInterpolation(WorkPackage &work_package) {
     this->pht->incrementReadCounter(roundKey(rounded_key));
 #endif
 
+    const int cell_id = static_cast<int>(work_package.input[wp_i][0]);
+
+    if (!to_calc_cache.contains(cell_id)) {
+      const std::vector<std::int32_t> &to_calc = dht_instance.getKeyElements();
+      std::vector<std::int32_t> keep_indices;
+
+      for (std::size_t i = 0; i < to_calc.size(); i++) {
+        if (!std::isnan(work_package.input[wp_i][to_calc[i]])) {
+          keep_indices.push_back(to_calc[i]);
+        }
+      }
+
+      to_calc_cache[cell_id] = keep_indices;
+    }
+
     double start_fc = MPI_Wtime();
 
     work_package.output[wp_i] =
-        f_interpolate(dht_instance.getKeyElements(), work_package.input[wp_i],
+        f_interpolate(to_calc_cache[cell_id], work_package.input[wp_i],
                       pht_result.in_values, pht_result.out_values);
 
     if (hooks.interp_post.isValid()) {
       NamedVector<double> nv_result(this->out_names, work_package.output[wp_i]);
-      if (hooks.interp_post(nv_result)) {
+      if (Rcpp::as<bool>(hooks.interp_post(nv_result))) {
         interp_result.status[wp_i] = INSUFFICIENT_DATA;
         continue;
       }

src/Chemistry/SurrogateModels/LookupKey.hpp

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

src/Chemistry/SurrogateModels/ProximityHashTable.cpp

@@ -1,4 +1,3 @@
-//  Time-stamp: "Last modified 2023-08-15 14:50:59 mluebke"
 #include "Interpolation.hpp"
 
 #include "DHT_Wrapper.hpp"

src/Chemistry/SurrogateModels/Rounding.hpp

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

src/Chemistry/WorkerFunctions.cpp

@@ -9,7 +9,6 @@
 #include <cstdint>
 #include <iomanip>
 #include <iostream>
-#include <map>
 #include <mpi.h>
 #include <stdexcept>
 #include <string>
@@ -47,6 +46,16 @@ void poet::ChemistryModule::WorkerLoop() {
     switch (func_type) {
     case CHEM_FIELD_INIT: {
       ChemBCast(&this->prop_count, 1, MPI_UINT32_T);
+      if (this->ai_surrogate_enabled) {
+        this->ai_surrogate_validity_vector.resize(
+            this->n_cells); // resize statt reserve?
+      }
+      break;
+    }
+    case CHEM_AI_BCAST_VALIDITY: {
+      // Receive the index vector of valid ai surrogate predictions
+      MPI_Bcast(&this->ai_surrogate_validity_vector.front(), this->n_cells,
+                MPI_INT, 0, this->group_comm);
       break;
     }
     case CHEM_WORK_LOOP: {
@@ -118,7 +127,7 @@ void poet::ChemistryModule::WorkerDoWork(MPI_Status &probe_status,
   uint32_t iteration;
   double dt;
   double current_sim_time;
-
+  uint32_t wp_start_index;
   int count = double_count;
   std::vector<double> mpi_buffer(count);
 
@@ -144,8 +153,8 @@ void poet::ChemistryModule::WorkerDoWork(MPI_Status &probe_status,
   // current simulation time ('age' of simulation)
   current_sim_time = mpi_buffer[count + 3];
 
-  /* 4th double value is currently a placeholder */
-  // placeholder = mpi_buffer[count+4];
+  // current work package start location in field
+  wp_start_index = mpi_buffer[count + 4];
 
   for (std::size_t wp_i = 0; wp_i < s_curr_wp.size; wp_i++) {
     s_curr_wp.input[wp_i] =
@@ -170,6 +179,15 @@ void poet::ChemistryModule::WorkerDoWork(MPI_Status &probe_status,
     interp->tryInterpolation(s_curr_wp);
   }
 
+  if (this->ai_surrogate_enabled) {
+    // Map valid predictions from the ai surrogate in the workpackage
+    for (int i = 0; i < s_curr_wp.size; i++) {
+      if (this->ai_surrogate_validity_vector[wp_start_index + i] == 1) {
+        s_curr_wp.mapping[i] = CHEM_AISURR;
+      }
+    }
+  }
+
   phreeqc_time_start = MPI_Wtime();
 
   WorkerRunWorkPackage(s_curr_wp, current_sim_time, dt);
@@ -228,6 +246,17 @@ void poet::ChemistryModule::WorkerPostIter(MPI_Status &prope_status,
           << std::setw(this->file_pad) << iteration << ".pht";
       interp->dumpPHTState(out.str());
     }
+
+    const auto max_mean_idx =
+        DHT_get_used_idx_factor(this->interp->getDHTObject(), 1);
+
+    if (max_mean_idx >= 2) {
+      DHT_flush(this->interp->getDHTObject());
+      DHT_flush(this->dht->getDHT());
+      if (this->comm_rank == 2) {
+        std::cout << "Flushed both DHT and PHT!\n\n";
+      }
+    }
   }
 
   RInsidePOET::getInstance().parseEvalQ("gc()");
@@ -282,28 +311,19 @@ void poet::ChemistryModule::WorkerRunWorkPackage(WorkPackage &work_package,
                                                  double dSimTime,
                                                  double dTimestep) {
 
+  std::vector<std::vector<double>> inout_chem = work_package.input;
+  std::vector<std::size_t> to_ignore;
+
   for (std::size_t wp_id = 0; wp_id < work_package.size; wp_id++) {
     if (work_package.mapping[wp_id] != CHEM_PQC) {
-      continue;
+      to_ignore.push_back(wp_id);
     }
+  }
+  this->pqc_runner->run(inout_chem, dTimestep, to_ignore);
 
-    auto curr_input = work_package.input[wp_id];
-    const auto pqc_id = static_cast<int>(curr_input[0]);
-
-    auto &phreeqc_instance = this->phreeqc_instances[pqc_id];
-    work_package.output[wp_id] = work_package.input[wp_id];
-
-    curr_input.erase(std::remove_if(curr_input.begin(), curr_input.end(),
-                                    [](double d) { return std::isnan(d); }),
-                     curr_input.end());
-
-    phreeqc_instance->runCell(curr_input, dTimestep);
-
-    std::size_t output_index = 0;
-    for (std::size_t i = 0; i < work_package.output[wp_id].size(); i++) {
-      if (!std::isnan(work_package.output[wp_id][i])) {
-        work_package.output[wp_id][i] = curr_input[output_index++];
-      }
+  for (std::size_t wp_id = 0; wp_id < work_package.size; wp_id++) {
+    if (work_package.mapping[wp_id] == CHEM_PQC) {
+      work_package.output[wp_id] = inout_chem[wp_id];
     }
   }
 }

src/Init/ChemistryInit.cpp

@@ -2,16 +2,14 @@
 
 #include <Rcpp.h>
 #include <cstddef>
+#include <fstream>
+#include <string>
+#include <utility>
 #include <vector>
 
 namespace poet {
 
 void InitialList::initChemistry(const Rcpp::List &chem) {
-  this->pqc_solutions = std::vector<std::string>(
-      this->transport_names.begin() + 3, this->transport_names.end());
-
-  this->pqc_solution_primaries = this->phreeqc->getSolutionPrimaries();
-
   if (chem.containsElementNamed("dht_species")) {
     this->dht_species = Rcpp::as<NamedVector<uint32_t>>(chem["dht_species"]);
   }
@@ -32,6 +30,30 @@ void InitialList::initChemistry(const Rcpp::List &chem) {
     }
   }
 
+  if (chem.containsElementNamed("ai_surrogate_input_script")) {
+    std::string ai_surrogate_input_script_path = chem["ai_surrogate_input_script"];
+    ai_surrogate_input_script_path = Rcpp::as<std::string>(Rcpp::Function("normalizePath")(Rcpp::wrap(ai_surrogate_input_script_path)));
+    
+    // Copying the entire script for the init file 
+    std::ifstream file(ai_surrogate_input_script_path);
+    if (!file.is_open()) {
+        // print error message and return
+        Rcpp::Rcerr << "AI surrogate input script was not found at: " << ai_surrogate_input_script_path << std::endl;
+    }
+
+    std::stringstream buffer;
+    buffer << file.rdbuf();
+    std::string fileContent = buffer.str();
+    file.close();
+    
+    // Get base path
+    ai_surrogate_input_script_path = ai_surrogate_input_script_path.substr(0, ai_surrogate_input_script_path.find_last_of('/') + 1);
+    // Add the filepath as a global variable in R to enable relative filepaths in the R script
+    fileContent += "\nai_surrogate_base_path <- \"" + ai_surrogate_input_script_path + "\"";
+    
+    this->ai_surrogate_input_script = fileContent;
+  }
+
   this->field_header =
       Rcpp::as<std::vector<std::string>>(this->initial_grid.names());
   this->field_header.erase(this->field_header.begin());
@@ -45,26 +67,22 @@ InitialList::ChemistryInit InitialList::getChemistryInit() const {
   // chem_init.field_header = this->field_header;
 
   chem_init.database = database;
+  chem_init.pqc_script = pqc_script;
+  chem_init.pqc_ids = pqc_ids;
+  chem_init.with_h0_o0 = with_h0_o0;
+  chem_init.with_redox = with_redox;
   // chem_init.pqc_scripts = pqc_scripts;
   // chem_init.pqc_ids = pqc_ids;
 
-  for (std::size_t i = 0; i < pqc_scripts.size(); i++) {
-    POETInitCell cell = {
-        pqc_solutions,
-        pqc_solution_primaries,
-        Rcpp::as<std::vector<std::string>>(pqc_exchanger[i]),
-        Rcpp::as<std::vector<std::string>>(pqc_kinetics[i]),
-        Rcpp::as<std::vector<std::string>>(pqc_equilibrium[i]),
-        Rcpp::as<std::vector<std::string>>(pqc_surface_comps[i]),
-        Rcpp::as<std::vector<std::string>>(pqc_surface_charges[i])};
-
-    chem_init.pqc_config[pqc_ids[i]] = {database, pqc_scripts[i], cell};
-  }
+  // for (std::size_t i = 0; i < pqc_ids.size(); ++i) {
+  //   chem_init.pqc_input[pqc_ids[i]] = pqc_scripts[i];
+  // }
 
   // chem_init.pqc_sol_order = pqc_solutions;
 
   chem_init.dht_species = dht_species;
   chem_init.interp_species = interp_species;
+  chem_init.ai_surrogate_input_script = ai_surrogate_input_script;
 
   if (this->chem_hooks.size() > 0) {
     if (this->chem_hooks.containsElementNamed("dht_fill")) {

src/Init/DiffusionInit.cpp

@@ -1,4 +1,3 @@
-#include <algorithm>
 #include <tug/Boundary.hpp>
 // leave above Rcpp includes, as eigen seem to have problems with a preceding
 // Rcpp include
@@ -8,7 +7,7 @@
 
 #include "DataStructures/Field.hpp"
 #include "InitialList.hpp"
-#include "PhreeqcInit.hpp"
+#include "PhreeqcMatrix.hpp"
 
 #include <Rcpp/Function.h>
 #include <Rcpp/proxy/ProtectedProxy.h>
@@ -27,9 +26,9 @@ namespace poet {
 enum SEXP_TYPE { SEXP_IS_LIST, SEXP_IS_VEC };
 
 const std::map<std::uint8_t, std::string> tug_side_mapping = {
-    {tug::BC_SIDE_RIGHT, "E"},
-    {tug::BC_SIDE_LEFT, "W"},
-    {tug::BC_SIDE_TOP, "N"},
+    {tug::BC_SIDE_RIGHT	, "E"},
+    {tug::BC_SIDE_LEFT	, "W"},
+    {tug::BC_SIDE_TOP	, "N"},
     {tug::BC_SIDE_BOTTOM, "S"}};
 
 static std::vector<TugType> colMajToRowMaj(const Rcpp::NumericVector &vec,
@@ -53,100 +52,104 @@ static std::vector<TugType> colMajToRowMaj(const Rcpp::NumericVector &vec,
   }
 }
 
-static std::vector<std::string>
-extend_transport_names(std::unique_ptr<PhreeqcInit> &phreeqc,
-                       const Rcpp::List &boundaries_list,
-                       const Rcpp::List &inner_boundaries,
-                       const std::vector<std::string> &old_trans_names) {
+// static std::vector<std::string>
+// extend_transport_names(std::unique_ptr<PhreeqcInit> &phreeqc,
+//                        const Rcpp::List &boundaries_list,
+//                        const Rcpp::List &inner_boundaries,
+//                        const std::vector<std::string> &old_trans_names) {
 
-  std::vector<std::string> transport_names = old_trans_names;
-  std::set<int> constant_pqc_ids;
+//   std::vector<std::string> transport_names = old_trans_names;
+//   std::set<int> constant_pqc_ids;
 
-  for (const auto &side : tug_side_mapping) {
-    if (!boundaries_list.containsElementNamed(side.second.c_str())) {
-      continue;
-    }
+//   for (const auto &side : tug_side_mapping) {
+//     if (!boundaries_list.containsElementNamed(side.second.c_str())) {
+//       continue;
+//     }
 
-    Rcpp::List mapping = boundaries_list[side.second];
+//     Rcpp::List mapping = boundaries_list[side.second];
 
-    const Rcpp::NumericVector cells = mapping["cell"];
-    const Rcpp::NumericVector values = mapping["sol_id"];
-    const Rcpp::CharacterVector type_str = mapping["type"];
+//     const Rcpp::NumericVector cells = mapping["cell"];
+//     const Rcpp::NumericVector values = mapping["sol_id"];
+//     const Rcpp::CharacterVector type_str = mapping["type"];
 
-    if (cells.size() != values.size()) {
-      throw std::runtime_error("Boundary [" + side.second +
-                               "] cells and values are not the same "
-                               "length");
-    }
+//     if (cells.size() != values.size()) {
+//       throw std::runtime_error("Boundary [" + side.second +
+//                                "] cells and values are not the same "
+//                                "length");
+//     }
 
-    for (auto i = 0; i < cells.size(); i++) {
-      if (type_str[i] == "constant") {
-        constant_pqc_ids.insert(static_cast<std::size_t>(values[i]));
-      }
-    }
-  }
+//     for (auto i = 0; i < cells.size(); i++) {
+//       if (type_str[i] == "constant") {
+//         constant_pqc_ids.insert(static_cast<std::size_t>(values[i]));
+//       }
+//     }
+//   }
 
-  if (inner_boundaries.size() > 0) {
-    const Rcpp::NumericVector values = inner_boundaries["sol_id"];
-    for (auto i = 0; i < values.size(); i++) {
-      constant_pqc_ids.insert(static_cast<std::size_t>(values[i]));
-    }
-  }
+//   if (inner_boundaries.size() > 0) {
+//     const Rcpp::NumericVector values = inner_boundaries["sol_id"];
+//     for (auto i = 0; i < values.size(); i++) {
+//       constant_pqc_ids.insert(static_cast<std::size_t>(values[i]));
+//     }
+//   }
 
-  if (!constant_pqc_ids.empty()) {
-    constexpr std::size_t keep_h_o_charge = 3;
+//   if (!constant_pqc_ids.empty()) {
+//     constexpr std::size_t keep_h_o_charge = 3;
 
-    for (const auto &pqc_id : constant_pqc_ids) {
-      const auto solution_names = phreeqc->getSolutionNames(pqc_id);
+//     for (const auto &pqc_id : constant_pqc_ids) {
+//       const auto solution_names = phreeqc->getSolutionNames(pqc_id);
 
-      // add those strings which are not already in the transport_names
-      for (const auto &name : solution_names) {
-        if (std::find(transport_names.begin(), transport_names.end(), name) ==
-            transport_names.end()) {
-          auto position =
-              std::lower_bound(transport_names.begin() + keep_h_o_charge,
-                               transport_names.end(), name);
+//       // add those strings which are not already in the transport_names
+//       for (const auto &name : solution_names) {
+//         if (std::find(transport_names.begin(), transport_names.end(), name)
+//         ==
+//             transport_names.end()) {
+//           auto position =
+//               std::lower_bound(transport_names.begin() + keep_h_o_charge,
+//                                transport_names.end(), name);
 
-          transport_names.insert(position, name);
-        }
-      }
-    }
-  }
+//           transport_names.insert(position, name);
+//         }
+//       }
+//     }
+//   }
 
-  return transport_names;
-}
+//   return transport_names;
+// }
 
-static Rcpp::List
-extend_initial_grid(const Rcpp::List &initial_grid,
-                    const std::vector<std::string> &transport_names) {
-  // std::vector<std::string> names_to_add(transport_names.begin() + old_size,
-  //                                       transport_names.end());
+// static Rcpp::List
+// extend_initial_grid(const Rcpp::List &initial_grid,
+//                     const std::vector<std::string> &transport_names) {
+//   // std::vector<std::string> names_to_add(transport_names.begin() +
+//   old_size,
+//   //                                       transport_names.end());
 
-  Rcpp::Function extend_grid_R("add_missing_transport_species");
+//   Rcpp::Function extend_grid_R("add_missing_transport_species");
 
-  return extend_grid_R(initial_grid, Rcpp::wrap(transport_names));
-}
+//   return extend_grid_R(initial_grid, Rcpp::wrap(transport_names));
+// }
 
 std::pair<Rcpp::List, Rcpp::List>
 InitialList::resolveBoundaries(const Rcpp::List &boundaries_list,
-                               const Rcpp::List &inner_boundaries) {
+                               const Rcpp::List &inner_boundaries,
+                               const PhreeqcMatrix &phreeqc_mat) {
   Rcpp::List bound_list;
   Rcpp::List inner_bound;
   Rcpp::Function resolve_R("resolve_pqc_bound");
 
-  const std::size_t old_transport_size = this->transport_names.size();
+  // const std::size_t old_transport_size = this->transport_names.size();
 
-  this->transport_names = extend_transport_names(
-      this->phreeqc, boundaries_list, inner_boundaries, this->transport_names);
+  // this->transport_names = extend_transport_names(
+  //     this->phreeqc, boundaries_list, inner_boundaries,
+  //     this->transport_names);
 
-  const std::size_t new_transport_size = this->transport_names.size();
+  // const std::size_t new_transport_size = this->transport_names.size();
 
-  if (old_transport_size != new_transport_size) {
-    this->initial_grid =
-        extend_initial_grid(this->initial_grid, this->transport_names);
-  }
+  // if (old_transport_size != new_transport_size) {
+  //   this->initial_grid =
+  //       extend_initial_grid(this->initial_grid, this->transport_names);
+  // }
 
-  for (const auto &species : this->transport_names) {
+  for (const auto &species_name : this->transport_names) {
     Rcpp::List spec_list;
 
     for (const auto &side : tug_side_mapping) {
@@ -163,8 +166,8 @@ InitialList::resolveBoundaries(const Rcpp::List &boundaries_list,
       if (boundaries_list.containsElementNamed(side.second.c_str())) {
         const Rcpp::List mapping = boundaries_list[side.second];
 
-        const Rcpp::NumericVector cells = mapping["cell"];
-        const Rcpp::NumericVector values = mapping["sol_id"];
+        const Rcpp::IntegerVector cells  = mapping["cell"];   // MDL 2024-06-13 
+        const Rcpp::IntegerVector values = mapping["sol_id"]; // MDL 
         const Rcpp::CharacterVector type_str = mapping["type"];
 
         if (cells.size() != values.size()) {
@@ -179,8 +182,11 @@ InitialList::resolveBoundaries(const Rcpp::List &boundaries_list,
             c_type[c_id] = tug::BC_TYPE_CLOSED;
           } else if (type_str[i] == "constant") {
             c_type[c_id] = tug::BC_TYPE_CONSTANT;
-            c_value[c_id] = Rcpp::as<TugType>(
-                resolve_R(this->phreeqc_mat, Rcpp::wrap(species), values[i]));
+            c_value[c_id] =
+                static_cast<TugType>(phreeqc_mat(values[i], species_name));
+            // c_value[c_id] = Rcpp::as<TugType>(
+            //     resolve_R(this->phreeqc_mat, Rcpp::wrap(species_name),
+            //     values[i]));
           } else {
             throw std::runtime_error("Unknown boundary type");
           }
@@ -191,7 +197,7 @@ InitialList::resolveBoundaries(const Rcpp::List &boundaries_list,
           Rcpp::Named("type") = c_type, Rcpp::Named("value") = c_value);
     }
 
-    bound_list[species] = spec_list;
+    bound_list[species_name] = spec_list;
 
     if (inner_boundaries.size() > 0) {
       const Rcpp::NumericVector &inner_row_vec = inner_boundaries["row"];
@@ -211,11 +217,14 @@ InitialList::resolveBoundaries(const Rcpp::List &boundaries_list,
       for (std::size_t i = 0; i < inner_row_vec.size(); i++) {
         rows.push_back(inner_row_vec[i] - 1);
         cols.push_back(inner_col_vec[i] - 1);
-        c_value.push_back(Rcpp::as<TugType>(resolve_R(
-            this->phreeqc_mat, Rcpp::wrap(species), inner_pqc_id_vec[i])));
+        c_value.push_back(static_cast<TugType>(
+            phreeqc_mat(inner_pqc_id_vec[i], species_name)));
+        // c_value.push_back(Rcpp::as<TugType>(resolve_R(
+        //     this->phreeqc_mat, Rcpp::wrap(species_name),
+        //     inner_pqc_id_vec[i])));
       }
 
-      inner_bound[species] =
+      inner_bound[species_name] =
           Rcpp::List::create(Rcpp::Named("row") = Rcpp::wrap(rows),
                              Rcpp::Named("col") = Rcpp::wrap(cols),
                              Rcpp::Named("value") = Rcpp::wrap(c_value));
@@ -276,7 +285,8 @@ static Rcpp::List parseAlphas(const SEXP &input,
 
   return out_list;
 }
-void InitialList::initDiffusion(const Rcpp::List &diffusion_input) {
+void InitialList::initDiffusion(const Rcpp::List &diffusion_input,
+                                const PhreeqcMatrix &phreeqc) {
   Rcpp::List boundaries;
   Rcpp::List inner_boundaries;
 
@@ -298,7 +308,7 @@ void InitialList::initDiffusion(const Rcpp::List &diffusion_input) {
       diffusion_input[DIFFU_MEMBER_STR(DiffusionMembers::ALPHA_Y)];
 
   const auto resolved_boundaries =
-      resolveBoundaries(boundaries, inner_boundaries);
+      resolveBoundaries(boundaries, inner_boundaries, phreeqc);
   this->boundaries = resolved_boundaries.first;
   this->inner_boundaries = resolved_boundaries.second;
 
@@ -403,4 +413,4 @@ InitialList::DiffusionInit InitialList::getDiffusionInit() const {
   return diff_init;
 }
 
-} // namespace poet
+} // namespace poet

src/Init/GridInit.cpp

@@ -1,13 +1,19 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+/*
+ * GridInit.cpp - Implementation of grid initialization
+ * Copyright (C) 2024-2025 Max Luebke (University of Potsdam)
+ */
+
 #include "InitialList.hpp"
 
-#include "PhreeqcInit.hpp"
+#include "PhreeqcMatrix.hpp"
 
 #include <RInside.h>
 #include <Rcpp/Function.h>
+#include <Rcpp/vector/Matrix.h>
 #include <Rcpp/vector/instantiation.h>
-#include <cstdint>
+#include <fstream>
 #include <map>
-#include <memory>
 #include <regex>
 #include <sstream>
 #include <string>
@@ -15,38 +21,6 @@
 
 namespace poet {
 
-static Rcpp::NumericMatrix
-pqcScriptToGrid(std::unique_ptr<PhreeqcInit> &phreeqc, RInside &R) {
-  PhreeqcInit::PhreeqcMat phreeqc_mat = phreeqc->getPhreeqcMat();
-
-  // add "id" to the front of the column names
-
-  const std::size_t ncols = phreeqc_mat.names.size();
-  const std::size_t nrows = phreeqc_mat.values.size();
-
-  phreeqc_mat.names.insert(phreeqc_mat.names.begin(), "ID");
-
-  Rcpp::NumericMatrix mat(nrows, ncols + 1);
-
-  for (std::size_t i = 0; i < nrows; i++) {
-    mat(i, 0) = phreeqc_mat.ids[i];
-    for (std::size_t j = 0; j < ncols; ++j) {
-      mat(i, j + 1) = phreeqc_mat.values[i][j];
-    }
-  }
-
-  Rcpp::colnames(mat) = Rcpp::wrap(phreeqc_mat.names);
-
-  return mat;
-}
-
-static inline Rcpp::List matToGrid(RInside &R, const Rcpp::NumericMatrix &mat,
-                                   const Rcpp::NumericMatrix &grid) {
-  Rcpp::Function pqc_to_grid_R("pqc_to_grid");
-
-  return pqc_to_grid_R(mat, grid);
-}
-
 static inline std::map<int, std::string>
 replaceRawKeywordIDs(std::map<int, std::string> raws) {
   for (auto &raw : raws) {
@@ -60,25 +34,6 @@ replaceRawKeywordIDs(std::map<int, std::string> raws) {
   return raws;
 }
 
-static inline uint32_t getSolutionCount(std::unique_ptr<PhreeqcInit> &phreeqc,
-                                        const Rcpp::List &initial_grid) {
-  PhreeqcInit::ModulesArray mod_array;
-  Rcpp::Function unique_R("unique");
-
-  std::vector<int> row_ids =
-      Rcpp::as<std::vector<int>>(unique_R(initial_grid["ID"]));
-
-  // std::vector<std::uint32_t> sizes_vec(sizes.begin(), sizes.end());
-
-  // Rcpp::Function modify_sizes("modify_module_sizes");
-  // sizes_vec = Rcpp::as<std::vector<std::uint32_t>>(
-  //     modify_sizes(sizes_vec, phreeqc_mat, initial_grid));
-
-  // std::copy(sizes_vec.begin(), sizes_vec.end(), sizes.begin());
-
-  return phreeqc->getModuleSizes(row_ids)[POET_SOL];
-}
-
 static std::string readFile(const std::string &path) {
   std::string string_rpath(PATH_MAX, '\0');
 
@@ -98,10 +53,26 @@ static std::string readFile(const std::string &path) {
   return buffer.str();
 }
 
-void InitialList::prepareGrid(const Rcpp::List &grid_input) {
-  // parse input values
-  Rcpp::Function unique_R("unique");
+static Rcpp::List expandGrid(const PhreeqcMatrix &pqc_mat,
+                             const std::vector<int> unique_ids,
+                             const Rcpp::IntegerMatrix &grid_def) {
 
+  PhreeqcMatrix subset_pqc_mat = pqc_mat.subset(unique_ids);
+
+  PhreeqcMatrix::STLExport phreeqc_mat =
+      subset_pqc_mat.get(PhreeqcMatrix::VectorExportType::COLUMN_MAJOR);
+
+  const std::size_t ncols = phreeqc_mat.names.size();
+  const std::size_t nrows = phreeqc_mat.values.size() / ncols;
+
+  Rcpp::NumericMatrix pqc_mat_R(nrows, ncols, phreeqc_mat.values.begin());
+  Rcpp::colnames(pqc_mat_R) = Rcpp::wrap(phreeqc_mat.names);
+
+  return Rcpp::Function("pqc_to_grid")(pqc_mat_R, grid_def);
+}
+
+PhreeqcMatrix InitialList::prepareGrid(const Rcpp::List &grid_input) {
+  // parse input values
   std::string script;
   std::string database;
 
@@ -126,8 +97,9 @@ void InitialList::prepareGrid(const Rcpp::List &grid_input) {
   }
 
   this->database = database;
+  this->pqc_script = script;
 
-  Rcpp::NumericMatrix grid_def =
+  Rcpp::IntegerMatrix grid_def =
       grid_input[GRID_MEMBER_STR(GridMembers::GRID_DEF)];
   Rcpp::NumericVector grid_size =
       grid_input[GRID_MEMBER_STR(GridMembers::GRID_SIZE)];
@@ -156,35 +128,77 @@ void InitialList::prepareGrid(const Rcpp::List &grid_input) {
     throw std::runtime_error("Grid size must be positive.");
   }
 
-  this->phreeqc = std::make_unique<PhreeqcInit>(database, script);
+  bool with_redox =
+      grid_input.containsElementNamed(
+          GRID_MEMBER_STR(GridMembers::PQC_WITH_REDOX))
+          ? Rcpp::as<bool>(
+                grid_input[GRID_MEMBER_STR(GridMembers::PQC_WITH_REDOX)])
+          : false;
 
-  this->phreeqc_mat = pqcScriptToGrid(phreeqc, R);
-  this->initial_grid = matToGrid(R, this->phreeqc_mat, grid_def);
+  bool with_h0_o0 =
+      grid_input.containsElementNamed(
+          GRID_MEMBER_STR(GridMembers::PQC_WITH_H0_O0))
+          ? Rcpp::as<bool>(
+                grid_input[GRID_MEMBER_STR(GridMembers::PQC_WITH_H0_O0)])
+          : false;
 
-  const uint32_t solution_count = getSolutionCount(phreeqc, this->initial_grid);
-
-  std::vector<std::string> colnames =
-      Rcpp::as<std::vector<std::string>>(this->initial_grid.names());
-
-  this->transport_names = std::vector<std::string>(
-      colnames.begin() + 1,
-      colnames.begin() + 1 + solution_count); // skip ID
-
-  std::map<int, std::string> pqc_raw_dumps;
-
-  pqc_raw_dumps = replaceRawKeywordIDs(phreeqc->raw_dumps());
-
-  this->pqc_ids =
-      Rcpp::as<std::vector<int>>(unique_R(this->initial_grid["ID"]));
-
-  for (const auto &id : this->pqc_ids) {
-    this->pqc_scripts.push_back(pqc_raw_dumps[id]);
-    this->pqc_exchanger.push_back(phreeqc->getExchanger(id));
-    this->pqc_kinetics.push_back(phreeqc->getKineticsNames(id));
-    this->pqc_equilibrium.push_back(phreeqc->getEquilibriumNames(id));
-    this->pqc_surface_comps.push_back(phreeqc->getSurfaceCompNames(id));
-    this->pqc_surface_charges.push_back(phreeqc->getSurfaceChargeNames(id));
+  if (with_h0_o0 && !with_redox) {
+    throw std::runtime_error(
+        "Output of H(0) and O(0) can only be used with redox.");
   }
+
+  this->with_h0_o0 = with_h0_o0;
+  this->with_redox = with_redox;
+
+  PhreeqcMatrix pqc_mat =
+      PhreeqcMatrix(database, script, with_h0_o0, with_redox);
+
+  this->transport_names = pqc_mat.getMatrixTransported();
+
+  Rcpp::Function unique_R("unique");
+  Rcpp::Function sort_R("sort");
+
+  std::vector<int> unique_ids = Rcpp::as<std::vector<int>>(
+      unique_R(Rcpp::IntegerVector(grid_def.begin(), grid_def.end())));
+
+  this->initial_grid = expandGrid(pqc_mat, unique_ids, grid_def);
+
+  const auto pqc_raw_dumps = replaceRawKeywordIDs(pqc_mat.getDumpStringsPQI());
+
+  for (const auto &id : unique_ids) {
+    this->pqc_ids.push_back(id);
+  }
+
+  return pqc_mat;
+  // this->phreeqc_mat = pqcScriptToGrid(phreeqc, R);
+  // this->initial_grid = matToGrid(R, this->phreeqc_mat, grid_def);
+
+  // const uint32_t solution_count = getSolutionCount(phreeqc,
+  // this->initial_grid);
+
+  // std::vector<std::string> colnames =
+  //     Rcpp::as<std::vector<std::string>>(this->initial_grid.names());
+
+  // this->transport_names = std::vector<std::string>(
+  //     colnames.begin() + 1,
+  //     colnames.begin() + 1 + solution_count); // skip ID
+
+  // std::map<int, std::string> pqc_raw_dumps;
+
+  // pqc_raw_dumps = replaceRawKeywordIDs(phreeqc->raw_dumps());
+
+  // this->pqc_ids =
+  //     Rcpp::as<std::vector<int>>(unique_R(this->initial_grid["ID"]));
+
+  // for (const auto &id : this->pqc_ids) {
+  //   this->pqc_scripts.push_back(pqc_raw_dumps[id]);
+  //   // this->pqc_exchanger.push_back(phreeqc->getExchanger(id));
+  //   // this->pqc_kinetics.push_back(phreeqc->getKineticsNames(id));
+  //   // this->pqc_equilibrium.push_back(phreeqc->getEquilibriumNames(id));
+  //   // this->pqc_surface_comps.push_back(phreeqc->getSurfaceCompNames(id));
+  //   //
+  //   this->pqc_surface_charges.push_back(phreeqc->getSurfaceChargeNames(id));
+  // }
 }
 
-} // namespace poet
+} // namespace poet

src/Init/InitialList.cpp

@@ -1,5 +1,6 @@
 #include "InitialList.hpp"
 #include "DataStructures/NamedVector.hpp"
+#include "PhreeqcMatrix.hpp"
 #include <Rcpp/internal/wrap.h>
 #include <Rcpp/iostream/Rstreambuf.h>
 #include <Rcpp/proxy/ProtectedProxy.h>
@@ -10,8 +11,8 @@
 
 namespace poet {
 void InitialList::initializeFromList(const Rcpp::List &setup) {
-  prepareGrid(setup[grid_key]);
-  initDiffusion(setup[diffusion_key]);
+  PhreeqcMatrix phreeqc = prepareGrid(setup[grid_key]);
+  initDiffusion(setup[diffusion_key], phreeqc);
   initChemistry(setup[chemistry_key]);
 }
 
@@ -53,26 +54,30 @@ void InitialList::importList(const Rcpp::List &setup, bool minimal) {
   }
   this->database =
       Rcpp::as<std::string>(setup[static_cast<int>(ExportList::CHEM_DATABASE)]);
+  this->pqc_script = Rcpp::as<std::string>(
+      setup[static_cast<int>(ExportList::CHEM_PQC_SCRIPT)]);
+  this->with_h0_o0 =
+      Rcpp::as<bool>(setup[static_cast<int>(ExportList::CHEM_PQC_WITH_H0_O0)]);
+  this->with_redox =
+      Rcpp::as<bool>(setup[static_cast<int>(ExportList::CHEM_PQC_WITH_REDOX)]);
   this->field_header = Rcpp::as<std::vector<std::string>>(
       setup[static_cast<int>(ExportList::CHEM_FIELD_HEADER)]);
-  this->pqc_scripts = Rcpp::as<std::vector<std::string>>(
-      setup[static_cast<int>(ExportList::CHEM_PQC_SCRIPTS)]);
   this->pqc_ids = Rcpp::as<std::vector<int>>(
       setup[static_cast<int>(ExportList::CHEM_PQC_IDS)]);
-  this->pqc_solutions = Rcpp::as<std::vector<std::string>>(
-      setup[static_cast<int>(ExportList::CHEM_PQC_SOLUTIONS)]);
-  this->pqc_solution_primaries = Rcpp::as<std::vector<std::string>>(
-      setup[static_cast<int>(ExportList::CHEM_PQC_SOLUTION_PRIMARY)]);
-  this->pqc_exchanger =
-      Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_EXCHANGER)]);
-  this->pqc_kinetics =
-      Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_KINETICS)]);
-  this->pqc_equilibrium =
-      Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_EQUILIBRIUM)]);
-  this->pqc_surface_comps =
-      Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_SURFACE_COMPS)]);
-  this->pqc_surface_charges =
-      Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_SURFACE_CHARGES)]);
+  //   this->pqc_solutions = Rcpp::as<std::vector<std::string>>(
+  //       setup[static_cast<int>(ExportList::CHEM_PQC_SOLUTIONS)]);
+  //   this->pqc_solution_primaries = Rcpp::as<std::vector<std::string>>(
+  //       setup[static_cast<int>(ExportList::CHEM_PQC_SOLUTION_PRIMARY)]);
+  //   this->pqc_exchanger =
+  //       Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_EXCHANGER)]);
+  //   this->pqc_kinetics =
+  //       Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_KINETICS)]);
+  //   this->pqc_equilibrium =
+  //       Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_EQUILIBRIUM)]);
+  //   this->pqc_surface_comps =
+  //       Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_SURFACE_COMPS)]);
+  //   this->pqc_surface_charges =
+  //       Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_SURFACE_CHARGES)]);
 
   this->dht_species = NamedVector<uint32_t>(
       setup[static_cast<int>(ExportList::CHEM_DHT_SPECIES)]);
@@ -82,6 +87,8 @@ void InitialList::importList(const Rcpp::List &setup, bool minimal) {
 
   this->chem_hooks =
       Rcpp::as<Rcpp::List>(setup[static_cast<int>(ExportList::CHEM_HOOKS)]);
+
+  this->ai_surrogate_input_script = Rcpp::as<std::string>(setup[static_cast<int>(ExportList::AI_SURROGATE_INPUT_SCRIPT)]);
 }
 
 Rcpp::List InitialList::exportList() {
@@ -106,29 +113,34 @@ Rcpp::List InitialList::exportList() {
   out[static_cast<int>(ExportList::DIFFU_ALPHA_Y)] = this->alpha_y;
 
   out[static_cast<int>(ExportList::CHEM_DATABASE)] = Rcpp::wrap(this->database);
+  out[static_cast<int>(ExportList::CHEM_PQC_SCRIPT)] =
+      Rcpp::wrap(this->pqc_script);
+  out[static_cast<int>(ExportList::CHEM_PQC_WITH_H0_O0)] =
+      Rcpp::wrap(this->with_h0_o0);
+  out[static_cast<int>(ExportList::CHEM_PQC_WITH_REDOX)] =
+      Rcpp::wrap(this->with_redox);
   out[static_cast<int>(ExportList::CHEM_FIELD_HEADER)] =
       Rcpp::wrap(this->field_header);
-  out[static_cast<int>(ExportList::CHEM_PQC_SCRIPTS)] =
-      Rcpp::wrap(this->pqc_scripts);
   out[static_cast<int>(ExportList::CHEM_PQC_IDS)] = Rcpp::wrap(this->pqc_ids);
-  out[static_cast<int>(ExportList::CHEM_PQC_SOLUTIONS)] =
-      Rcpp::wrap(this->pqc_solutions);
-  out[static_cast<int>(ExportList::CHEM_PQC_SOLUTION_PRIMARY)] =
-      Rcpp::wrap(this->pqc_solution_primaries);
-  out[static_cast<int>(ExportList::CHEM_PQC_EXCHANGER)] =
-      Rcpp::wrap(this->pqc_exchanger);
-  out[static_cast<int>(ExportList::CHEM_PQC_KINETICS)] =
-      Rcpp::wrap(this->pqc_kinetics);
-  out[static_cast<int>(ExportList::CHEM_PQC_EQUILIBRIUM)] =
-      Rcpp::wrap(this->pqc_equilibrium);
-  out[static_cast<int>(ExportList::CHEM_PQC_SURFACE_COMPS)] =
-      Rcpp::wrap(this->pqc_surface_comps);
-  out[static_cast<int>(ExportList::CHEM_PQC_SURFACE_CHARGES)] =
-      Rcpp::wrap(this->pqc_surface_charges);
+  //   out[static_cast<int>(ExportList::CHEM_PQC_SOLUTIONS)] =
+  //       Rcpp::wrap(this->pqc_solutions);
+  //   out[static_cast<int>(ExportList::CHEM_PQC_SOLUTION_PRIMARY)] =
+  //       Rcpp::wrap(this->pqc_solution_primaries);
+  //   out[static_cast<int>(ExportList::CHEM_PQC_EXCHANGER)] =
+  //       Rcpp::wrap(this->pqc_exchanger);
+  //   out[static_cast<int>(ExportList::CHEM_PQC_KINETICS)] =
+  //       Rcpp::wrap(this->pqc_kinetics);
+  //   out[static_cast<int>(ExportList::CHEM_PQC_EQUILIBRIUM)] =
+  //       Rcpp::wrap(this->pqc_equilibrium);
+  //   out[static_cast<int>(ExportList::CHEM_PQC_SURFACE_COMPS)] =
+  //       Rcpp::wrap(this->pqc_surface_comps);
+  //   out[static_cast<int>(ExportList::CHEM_PQC_SURFACE_CHARGES)] =
+  //       Rcpp::wrap(this->pqc_surface_charges);
   out[static_cast<int>(ExportList::CHEM_DHT_SPECIES)] = this->dht_species;
   out[static_cast<int>(ExportList::CHEM_INTERP_SPECIES)] =
       Rcpp::wrap(this->interp_species);
   out[static_cast<int>(ExportList::CHEM_HOOKS)] = this->chem_hooks;
+  out[static_cast<int>(ExportList::AI_SURROGATE_INPUT_SCRIPT)] = this->ai_surrogate_input_script;
 
   return out;
 }

src/Init/InitialList.hpp

@@ -2,7 +2,7 @@
 
 #include "Base/RInsidePOET.hpp"
 #include "DataStructures/NamedVector.hpp"
-#include "POETInit.hpp"
+#include "PhreeqcMatrix.hpp"
 #include <Rcpp/vector/instantiation.h>
 #include <set>
 #include <tug/Boundary.hpp>
@@ -13,9 +13,7 @@
 #include <cstddef>
 #include <cstdint>
 
-#include <PhreeqcInit.hpp>
 #include <map>
-#include <memory>
 #include <string>
 #include <utility>
 #include <vector>
@@ -53,20 +51,16 @@ private:
     DIFFU_ALPHA_X,
     DIFFU_ALPHA_Y,
     CHEM_DATABASE,
-    CHEM_FIELD_HEADER,
-    CHEM_PQC_SCRIPTS,
+    CHEM_PQC_SCRIPT,
+    CHEM_PQC_WITH_H0_O0,
+    CHEM_PQC_WITH_REDOX,
     CHEM_PQC_IDS,
-    CHEM_PQC_SOLUTIONS,
-    CHEM_PQC_SOLUTION_PRIMARY,
-    CHEM_PQC_EXCHANGER,
-    CHEM_PQC_KINETICS,
-    CHEM_PQC_EQUILIBRIUM,
-    CHEM_PQC_SURFACE_COMPS,
-    CHEM_PQC_SURFACE_CHARGES,
+    CHEM_FIELD_HEADER,
     CHEM_DHT_SPECIES,
     CHEM_INTERP_SPECIES,
     CHEM_HOOKS,
-    ENUM_SIZE
+    AI_SURROGATE_INPUT_SCRIPT,
+    ENUM_SIZE // Hack: Last element of the enum to show enum size
   };
 
   // Grid members
@@ -75,6 +69,8 @@ private:
   enum class GridMembers {
     PQC_SCRIPT_STRING,
     PQC_SCRIPT_FILE,
+    PQC_WITH_REDOX,
+    PQC_WITH_H0_O0,
     PQC_DB_STRING,
     PQC_DB_FILE,
     GRID_DEF,
@@ -88,17 +84,18 @@ private:
       static_cast<std::size_t>(InitialList::GridMembers::ENUM_SIZE);
 
   static constexpr std::array<const char *, size_GridMembers>
-      GridMembersString = {"pqc_in_string",  "pqc_in_file", "pqc_db_string",
-                           "pqc_db_file",    "grid_def",    "grid_size",
-                           "constant_cells", "porosity"};
+      GridMembersString = {"pqc_in_string", "pqc_in_file",   "pqc_with_redox",
+                           "pqc_wth_h0_o0", "pqc_db_string", "pqc_db_file",
+                           "grid_def",      "grid_size",     "constant_cells",
+                           "porosity"};
 
   constexpr const char *GRID_MEMBER_STR(GridMembers member) const {
     return GridMembersString[static_cast<std::size_t>(member)];
   }
 
-  std::unique_ptr<PhreeqcInit> phreeqc;
+  // std::unique_ptr<PhreeqcMatrix> pqc_mat;
 
-  void prepareGrid(const Rcpp::List &grid_input);
+  PhreeqcMatrix prepareGrid(const Rcpp::List &grid_input);
 
   std::uint8_t dim{0};
 
@@ -113,9 +110,6 @@ private:
 
   Rcpp::List initial_grid;
 
-  // No export
-  Rcpp::NumericMatrix phreeqc_mat;
-
 public:
   struct DiffusionInit {
     using BoundaryMap =
@@ -168,10 +162,12 @@ private:
     return DiffusionMembersString[static_cast<std::size_t>(member)];
   }
 
-  void initDiffusion(const Rcpp::List &diffusion_input);
+  void initDiffusion(const Rcpp::List &diffusion_input,
+                     const PhreeqcMatrix &phreeqc);
   std::pair<Rcpp::List, Rcpp::List>
   resolveBoundaries(const Rcpp::List &boundaries_list,
-                    const Rcpp::List &inner_boundaries);
+                    const Rcpp::List &inner_boundaries,
+                    const PhreeqcMatrix &phreeqc_mat);
 
   Rcpp::List boundaries;
   Rcpp::List inner_boundaries;
@@ -189,21 +185,19 @@ private:
   std::vector<std::string> field_header;
 
   std::string database;
-  std::vector<std::string> pqc_scripts;
+  std::string pqc_script;
+  bool with_h0_o0{false};
+  bool with_redox{false};
+  // std::vector<std::string> pqc_scripts;
   std::vector<int> pqc_ids;
 
-  std::vector<std::string> pqc_solutions;
-  std::vector<std::string> pqc_solution_primaries;
-  Rcpp::List pqc_exchanger;
-  Rcpp::List pqc_kinetics;
-  Rcpp::List pqc_equilibrium;
-  Rcpp::List pqc_surface_comps;
-  Rcpp::List pqc_surface_charges;
-
   NamedVector<std::uint32_t> dht_species;
 
   NamedVector<std::uint32_t> interp_species;
 
+  // Path to R script that the user defines in the input file
+  std::string ai_surrogate_input_script;
+
   Rcpp::List chem_hooks;
 
   const std::set<std::string> hook_name_list{"dht_fill", "dht_fuzz",
@@ -211,10 +205,10 @@ private:
 
 public:
   struct ChemistryHookFunctions {
-    RHookFunction<bool> dht_fill;
-    RHookFunction<std::vector<double>> dht_fuzz;
-    RHookFunction<std::vector<std::size_t>> interp_pre;
-    RHookFunction<bool> interp_post;
+    poet::DEFunc dht_fill;
+    poet::DEFunc dht_fuzz;
+    poet::DEFunc interp_pre;
+    poet::DEFunc interp_post;
   };
 
   struct ChemistryInit {
@@ -223,16 +217,20 @@ public:
     // std::vector<std::string> field_header;
 
     std::string database;
-    // std::vector<std::string> pqc_scripts;
-    // std::vector<int> pqc_ids;
+    std::string pqc_script;
+    bool with_h0_o0;
+    bool with_redox;
+    std::vector<int> pqc_ids;
 
-    std::map<int, POETConfig> pqc_config;
+    // std::map<int, std::string> pqc_input;
 
     // std::vector<std::string> pqc_sol_order;
 
     NamedVector<std::uint32_t> dht_species;
     NamedVector<std::uint32_t> interp_species;
     ChemistryHookFunctions hooks;
+
+    std::string ai_surrogate_input_script;
   };
 
   ChemistryInit getChemistryInit() const;

src/Transport/DiffusionModule.cpp

@@ -1,97 +1,48 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
 /*
-** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
-** Potsdam)
-**
-** Copyright (C) 2018-2022 Marco De Lucia, Max Luebke (GFZ Potsdam)
-**
-** Copyright (C) 2023-2024 Marco De Lucia (GFZ Potsdam), Max Luebke (University
-** of Potsdam)
-**
-** POET is free software; you can redistribute it and/or modify it under the
-** terms of the GNU General Public License as published by the Free Software
-** Foundation; either version 2 of the License, or (at your option) any later
-** 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.
-*/
+ * DiffusionModule.cpp - Coupling module betweeen POET and tug
+ * Copyright (C) 2018-2025 Max Luebke (University of Potsdam)
+ */
 
 #include "DiffusionModule.hpp"
 
 #include "Base/Macros.hpp"
 #include "Init/InitialList.hpp"
 
-#include <Eigen/src/Core/Map.h>
-#include <Eigen/src/Core/Matrix.h>
-#include <Eigen/src/Core/util/Constants.h>
 #include <chrono>
-
-#include "tug/Boundary.hpp"
-#include "tug/Grid.hpp"
-#include "tug/Simulation.hpp"
-
 #include <string>
-#include <vector>
+#include <tug/Boundary.hpp>
+#include <tug/Diffusion.hpp>
 
 using namespace poet;
 
-static inline std::vector<TugType>
-MatrixToVec(const Eigen::MatrixX<TugType> &mat) {
-  std::vector<TugType> vec(mat.rows() * mat.cols());
-
-  for (std::uint32_t i = 0; i < mat.cols(); i++) {
-    for (std::uint32_t j = 0; j < mat.rows(); j++) {
-      vec[j * mat.cols() + i] = mat(j, i);
-    }
-  }
-
-  return vec;
-}
-
-static inline Eigen::MatrixX<TugType>
-VecToMatrix(const std::vector<TugType> &vec, std::uint32_t n_rows,
-            std::uint32_t n_cols) {
-  Eigen::MatrixX<TugType> mat(n_rows, n_cols);
-
-  for (std::uint32_t i = 0; i < n_cols; i++) {
-    for (std::uint32_t j = 0; j < n_rows; j++) {
-      mat(j, i) = vec[j * n_cols + i];
-    }
-  }
-
-  return mat;
-}
-
-// static constexpr double ZERO_MULTIPLICATOR = 10E-14;
-
 void DiffusionModule::simulate(double requested_dt) {
-  MSG("Starting diffusion ...");
+  // MSG("Starting diffusion ...");
   const auto start_diffusion_t = std::chrono::high_resolution_clock::now();
 
   const auto &n_rows = this->param_list.n_rows;
   const auto &n_cols = this->param_list.n_cols;
 
-  tug::Grid<TugType> grid(param_list.n_rows, param_list.n_cols);
-  tug::Boundary boundary(grid);
-
-  grid.setDomain(param_list.s_rows, param_list.s_cols);
-
-  tug::Simulation<TugType> sim(grid, boundary);
-  sim.setIterations(1);
-
   for (const auto &sol_name : this->param_list.transport_names) {
-    auto &species_conc = this->transport_field[sol_name];
+#if defined(POET_TUG_BTCS)
+    tug::Diffusion<TugType> diffusion_solver(
+        this->transport_field[sol_name].data(), n_rows, n_cols);
+#elif defined(POET_TUG_FTCS)
+    tug::Diffusion<TugType, tug::FTCS_APPROACH> diffusion_solver(
+        this->transport_field[sol_name].data(), n_rows, n_cols);
+#else
+#error "No valid approach defined"
+#endif
 
-    Eigen::MatrixX<TugType> conc = VecToMatrix(species_conc, n_rows, n_cols);
-    Eigen::MatrixX<TugType> alpha_x =
-        VecToMatrix(this->param_list.alpha_x[sol_name], n_rows, n_cols);
-    Eigen::MatrixX<TugType> alpha_y =
-        VecToMatrix(this->param_list.alpha_y[sol_name], n_rows, n_cols);
+    tug::RowMajMatMap<TugType> alpha_x_map(
+        this->param_list.alpha_x[sol_name].data(), n_rows, n_cols);
+    tug::RowMajMatMap<TugType> alpha_y_map(
+        this->param_list.alpha_y[sol_name].data(), n_rows, n_cols);
+
+    diffusion_solver.setAlphaX(alpha_x_map);
+    diffusion_solver.setAlphaY(alpha_y_map);
+
+    auto &boundary = diffusion_solver.getBoundaryConditions();
 
     boundary.deserialize(this->param_list.boundaries[sol_name]);
 
@@ -99,14 +50,8 @@ void DiffusionModule::simulate(double requested_dt) {
       boundary.setInnerBoundaries(this->param_list.inner_boundaries[sol_name]);
     }
 
-    grid.setAlpha(alpha_x, alpha_y);
-    grid.setConcentrations(conc);
-
-    sim.setTimestep(requested_dt);
-
-    sim.run();
-
-    species_conc = MatrixToVec(grid.getConcentrations());
+    diffusion_solver.setTimestep(requested_dt);
+    diffusion_solver.run();
   }
 
   const auto end_diffusion_t = std::chrono::high_resolution_clock::now();

src/initializer.cpp

@@ -1,7 +1,7 @@
 #include "Init/InitialList.hpp"
 #include "poet.hpp"
 
-#include "Base/argh.hpp"
+#include <CLI/CLI.hpp>
 
 #include <cstdlib>
 
@@ -11,31 +11,41 @@
 #include <string>
 
 int main(int argc, char **argv) {
-
-  argh::parser cmdl({"-o", "--output"});
-
-  cmdl.parse(argc, argv);
-
-  if (cmdl[{"-h", "--help"}] || cmdl.pos_args().size() != 2) {
-    std::cout << "Usage: " << argv[0]
-              << " [-o, --output output_file]"
-                 " [-s, --setwd] "
-                 " <script.R>"
-              << std::endl;
-    return EXIT_SUCCESS;
-  }
-
+  // initialize RIinside
   RInside R(argc, argv);
 
   R.parseEvalQ(init_r_library);
+  R.parseEvalQ(kin_r_library);
 
-  std::string input_script = cmdl.pos_args()[1];
+  // parse command line arguments
+  CLI::App app{"POET Initializer - Poet R script to POET qs/rds converter"};
+
+  std::string input_script;
+  app.add_option("PoetScript.R", input_script, "POET R script to convert")
+      ->required();
+
+  std::string output_file;
+  app.add_option("-n, --name", output_file,
+                 "Name of the output file without file extension");
+
+  bool setwd;
+  app.add_flag("-s, --setwd", setwd,
+               "Set working directory to the directory of the directory of the "
+               "input script")
+      ->default_val(false);
+
+  bool asRDS{false};
+  app.add_flag("-r, --rds", asRDS, "Save output as .rds")->default_val(false);
+
+  bool asQS{false};
+  app.add_flag("-q, --qs", asQS, "Save output as .qs")->default_val(false);
+
+  CLI11_PARSE(app, argc, argv);
+
+  // source the input script
   std::string normalized_path_script;
   std::string in_file_name;
 
-  std::string curr_path =
-      Rcpp::as<std::string>(Rcpp::Function("normalizePath")(Rcpp::wrap(".")));
-
   try {
     normalized_path_script =
         Rcpp::as<std::string>(Rcpp::Function("normalizePath")(input_script));
@@ -51,14 +61,26 @@ int main(int argc, char **argv) {
     return EXIT_FAILURE;
   }
 
-  std::string output_file;
+  // if output file is not specified, use the input file name
+  if (output_file.empty()) {
+    std::string curr_path =
+        Rcpp::as<std::string>(Rcpp::Function("normalizePath")(Rcpp::wrap(".")));
 
-  cmdl({"-o", "--output"},
-       curr_path + "/" +
-           in_file_name.substr(0, in_file_name.find_last_of('.')) + ".rds") >>
-      output_file;
+    output_file = curr_path + "/" +
+                  in_file_name.substr(0, in_file_name.find_last_of('.'));
+  }
 
-  if (cmdl[{"-s", "--setwd"}]) {
+  // append the correct file extension
+  if (asRDS) {
+    output_file += ".rds";
+  } else if (asQS) {
+    output_file += ".qs";
+  } else {
+    output_file += ".qs2";
+  }
+
+  // set working directory to the directory of the input script
+  if (setwd) {
     const std::string dir_path = Rcpp::as<std::string>(
         Rcpp::Function("dirname")(normalized_path_script));
 
@@ -71,13 +93,12 @@ int main(int argc, char **argv) {
 
   init.initializeFromList(setup);
 
-  // replace file extension by .rds
-  Rcpp::Function save("saveRDS");
-
-  save(init.exportList(), Rcpp::wrap(output_file));
+  // use the generic handler defined in kin_r_library.R
+  Rcpp::Function SaveRObj_R("SaveRObj");
+  SaveRObj_R(init.exportList(), Rcpp::wrap(output_file));
 
   std::cout << "Saved result to " << output_file << std::endl;
 
   //   parseGrid(R, grid, results);
   return EXIT_SUCCESS;
-}
+}

src/poet.cpp

@@ -4,7 +4,8 @@
 **
 ** Copyright (C) 2018-2022 Marco De Lucia, Max Luebke (GFZ Potsdam)
 **
-** Copyright (C) 2023-2024 Max Luebke (University of Potsdam)
+** Copyright (C) 2023-2024 Marco De Lucia (GFZ Potsdam), Max Luebke (University
+** of Potsdam)
 **
 ** POET is free software; you can redistribute it and/or modify it under the
 ** terms of the GNU General Public License as published by the Free Software
@@ -22,6 +23,7 @@
 
 #include "Base/Macros.hpp"
 #include "Base/RInsidePOET.hpp"
+#include "CLI/CLI.hpp"
 #include "Chemistry/ChemistryModule.hpp"
 #include "DataStructures/Field.hpp"
 #include "Init/InitialList.hpp"
@@ -32,17 +34,21 @@
 #include <Rcpp/DataFrame.h>
 #include <Rcpp/Function.h>
 #include <Rcpp/vector/instantiation.h>
+#include <algorithm>
+#include <array>
+#include <cstdint>
 #include <cstdlib>
 #include <memory>
 #include <mpi.h>
-#include <optional>
+#include <set>
 #include <string>
 
-#include "Base/argh.hpp"
+#include <CLI/CLI.hpp>
 
 #include <poet.hpp>
+#include <vector>
 
-using namespace std;
+// using namespace std;
 using namespace poet;
 using namespace Rcpp;
 
@@ -50,17 +56,23 @@ static int MY_RANK = 0;
 
 static std::unique_ptr<Rcpp::List> global_rt_setup;
 
-// we need some layz evaluation, as we can't define the functions before the R
-// runtime is initialized
-static std::optional<Rcpp::Function> master_init_R;
-static std::optional<Rcpp::Function> master_iteration_end_R;
-static std::optional<Rcpp::Function> store_setup_R;
+// we need some lazy evaluation, as we can't define the functions
+// before the R runtime is initialized
+static poet::DEFunc master_init_R;
+static poet::DEFunc master_iteration_end_R;
+// MDL: unused -> static poet::DEFunc store_setup_R;
+static poet::DEFunc ReadRObj_R;
+static poet::DEFunc SaveRObj_R;
+static poet::DEFunc source_R;
 
 static void init_global_functions(RInside &R) {
   R.parseEval(kin_r_library);
-  master_init_R = Rcpp::Function("master_init");
-  master_iteration_end_R = Rcpp::Function("master_iteration_end");
-  store_setup_R = Rcpp::Function("StoreSetup");
+  master_init_R = DEFunc("master_init");
+  master_iteration_end_R = DEFunc("master_iteration_end");
+  // MDL: unused -> store_setup_R = DEFunc("StoreSetup");
+  source_R = DEFunc("source");
+  ReadRObj_R = DEFunc("ReadRObj");
+  SaveRObj_R = DEFunc("SaveRObj");
 }
 
 // HACK: this is a step back as the order and also the count of fields is
@@ -78,78 +90,108 @@ static void init_global_functions(RInside &R) {
 
 enum ParseRet { PARSER_OK, PARSER_ERROR, PARSER_HELP };
 
-ParseRet parseInitValues(char **argv, RuntimeParameters &params) {
-  // initialize argh object
-  argh::parser cmdl(argv);
+int parseInitValues(int argc, char **argv, RuntimeParameters &params) {
 
-  // if user asked for help
-  if (cmdl[{"help", "h"}]) {
-    if (MY_RANK == 0) {
-      MSG("Todo");
-      MSG("See README.md for further information.");
-    }
+  CLI::App app{"POET - Potsdam rEactive Transport simulator"};
 
-    return ParseRet::PARSER_HELP;
-  }
-  // if positional arguments are missing
-  if (!cmdl(3)) {
-    if (MY_RANK == 0) {
-      ERRMSG("POET needs 3 positional arguments: ");
-      ERRMSG("1) the R script defining your simulation runtime.");
-      ERRMSG("2) the initial .rds file generated by poet_init.");
-      ERRMSG("3) the directory prefix where to save results and profiling");
-    }
-    return ParseRet::PARSER_ERROR;
-  }
-
-  // parse flags and check for unknown
-  for (const auto &option : cmdl.flags()) {
-    if (!(flaglist.find(option) != flaglist.end())) {
-      if (MY_RANK == 0) {
-        ERRMSG("Unrecognized option: " + option);
-        ERRMSG("Make sure to use available options. Exiting!");
-      }
-      return ParseRet::PARSER_ERROR;
-    }
-  }
-
-  // parse parameters and check for unknown
-  for (const auto &option : cmdl.params()) {
-    if (!(paramlist.find(option.first) != paramlist.end())) {
-      if (MY_RANK == 0) {
-        ERRMSG("Unrecognized option: " + option.first);
-        ERRMSG("Make sure to use available options. Exiting!");
-      }
-      return ParseRet::PARSER_ERROR;
-    }
-  }
-
-  params.print_progressbar = cmdl[{"P", "progress"}];
+  app.add_flag("-P,--progress", params.print_progress,
+               "Print progress bar during chemical simulation");
 
   /*Parse work package size*/
-  cmdl("work-package-size", CHEM_DEFAULT_WORK_PACKAGE_SIZE) >>
-      params.work_package_size;
+  app.add_option(
+         "-w,--work-package-size", params.work_package_size,
+         "Work package size to distribute to each worker for chemistry module")
+      ->check(CLI::PositiveNumber)
+      ->default_val(RuntimeParameters::WORK_PACKAGE_SIZE_DEFAULT);
 
   /* Parse DHT arguments */
-  params.use_dht = cmdl["dht"];
-  params.use_interp = cmdl["interp"];
+  auto *dht_group = app.add_option_group("DHT", "DHT related options");
+
+  dht_group->add_flag("--dht", params.use_dht, "Enable DHT");
+
   // cout << "CPP: DHT is " << ( dht_enabled ? "ON" : "OFF" ) << '\n';
 
-  cmdl("dht-size", CHEM_DHT_SIZE_PER_PROCESS_MB) >> params.dht_size;
+  dht_group
+      ->add_option("--dht-size", params.dht_size,
+                   "DHT size per process in Megabyte")
+      ->check(CLI::PositiveNumber)
+      ->default_val(RuntimeParameters::DHT_SIZE_DEFAULT);
   // cout << "CPP: DHT size per process (Byte) = " << dht_size_per_process <<
   // endl;
 
-  cmdl("dht-snaps", 0) >> params.dht_snaps;
+  dht_group->add_option(
+      "--dht-snaps", params.dht_snaps,
+      "Save snapshots of DHT to disk: \n0 = disabled (default)\n1 = After "
+      "simulation\n2 = After each iteration");
 
-  params.use_interp = cmdl["interp"];
-  cmdl("interp-size", 100) >> params.interp_size;
-  cmdl("interp-min", 5) >> params.interp_min_entries;
-  cmdl("interp-bucket-entries", 20) >> params.interp_bucket_entries;
+  auto *interp_group =
+      app.add_option_group("Interpolation", "Interpolation related options");
+
+  interp_group->add_flag("--interp", params.use_interp, "Enable interpolation");
+  interp_group
+      ->add_option("--interp-size", params.interp_size,
+                   "Size of the interpolation table in Megabyte")
+      ->check(CLI::PositiveNumber)
+      ->default_val(RuntimeParameters::INTERP_SIZE_DEFAULT);
+  interp_group
+      ->add_option("--interp-min", params.interp_min_entries,
+                   "Minimum number of entries in the interpolation table")
+      ->check(CLI::PositiveNumber)
+      ->default_val(RuntimeParameters::INTERP_MIN_ENTRIES_DEFAULT);
+  interp_group
+      ->add_option(
+          "--interp-bucket-entries", params.interp_bucket_entries,
+          "Maximum number of entries in each bucket of the interpolation table")
+      ->check(CLI::PositiveNumber)
+      ->default_val(RuntimeParameters::INTERP_BUCKET_ENTRIES_DEFAULT);
+
+  app.add_flag("--ai-surrogate", params.use_ai_surrogate,
+               "Enable AI surrogate for chemistry module");
+
+  app.add_flag("--rds", params.as_rds,
+               "Save output as .rds file instead of default .qs2");
+
+  app.add_flag("--qs", params.as_qs,
+               "Save output as .qs file instead of default .qs2");
+
+  std::string init_file;
+  std::string runtime_file;
+
+  app.add_option("runtime_file", runtime_file,
+                 "Runtime R script defining the simulation")
+      ->required()
+      ->check(CLI::ExistingFile);
+
+  app.add_option(
+         "init_file", init_file,
+         "Initial R script defining the simulation, produced by poet_init")
+      ->required()
+      ->check(CLI::ExistingFile);
+
+  app.add_option("out_dir", params.out_dir,
+                 "Output directory of the simulation")
+      ->required();
+
+  try {
+    app.parse(argc, argv);
+  } catch (const CLI::ParseError &e) {
+    app.exit(e);
+    return -1;
+  }
+
+  // set the output extension
+  params.out_ext = "qs2";
+  if (params.as_rds)
+    params.out_ext = "rds";
+  if (params.as_qs)
+    params.out_ext = "qs";
 
   if (MY_RANK == 0) {
     // MSG("Complete results storage is " + BOOL_PRINT(simparams.store_result));
+    MSG("Output format/extension is " + params.out_ext);
     MSG("Work Package Size: " + std::to_string(params.work_package_size));
     MSG("DHT is " + BOOL_PRINT(params.use_dht));
+    MSG("AI Surrogate is " + BOOL_PRINT(params.use_ai_surrogate));
 
     if (params.use_dht) {
       // MSG("DHT strategy is " + std::to_string(simparams.dht_strategy));
@@ -173,14 +215,6 @@ ParseRet parseInitValues(char **argv, RuntimeParameters &params) {
           std::to_string(params.interp_bucket_entries));
     }
   }
-
-  std::string init_file;
-  std::string runtime_file;
-
-  cmdl(1) >> runtime_file;
-  cmdl(2) >> init_file;
-  cmdl(3) >> params.out_dir;
-
   // chem_params.dht_outdir = out_dir;
 
   /* distribute information to R runtime */
@@ -202,16 +236,17 @@ ParseRet parseInitValues(char **argv, RuntimeParameters &params) {
   // R["dht_log"] = simparams.dht_log;
 
   try {
-    Rcpp::Function source("source");
-    Rcpp::Function readRDS("readRDS");
 
-    Rcpp::List init_params_ = readRDS(init_file);
+    Rcpp::List init_params_(ReadRObj_R(init_file));
     params.init_params = init_params_;
 
     global_rt_setup = std::make_unique<Rcpp::List>();
-    *global_rt_setup = source(runtime_file, Rcpp::Named("local", true));
+    *global_rt_setup = source_R(runtime_file, Rcpp::Named("local", true));
     *global_rt_setup = global_rt_setup->operator[]("value");
 
+    // MDL add "out_ext" for output format to R setup
+    (*global_rt_setup)["out_ext"] = params.out_ext;
+
     params.timesteps =
         Rcpp::as<std::vector<double>>(global_rt_setup->operator[]("timesteps"));
 
@@ -250,35 +285,108 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, const RuntimeParameters &params,
    * master for the following loop) */
   uint32_t maxiter = params.timesteps.size();
 
-  if (params.print_progressbar) {
+  if (params.print_progress) {
     chem.setProgressBarPrintout(true);
   }
+  R["TMP_PROPS"] = Rcpp::wrap(chem.getField().GetProps());
 
   /* SIMULATION LOOP */
-
   double dSimTime{0};
   for (uint32_t iter = 1; iter < maxiter + 1; iter++) {
     double start_t = MPI_Wtime();
 
     const double &dt = params.timesteps[iter - 1];
 
-    //  cout << "CPP: Next time step is " << dt << "[s]" << endl;
-    MSG("Next time step is " + std::to_string(dt) + " [s]");
+    std::cout << std::endl;
 
     /* displaying iteration number, with C++ and R iterator */
-    MSG("Going through iteration " + std::to_string(iter));
+    MSG("Going through iteration " + std::to_string(iter) + "/" +
+        std::to_string(maxiter));
+
+    MSG("Current time step is " + std::to_string(dt));
 
     /* run transport */
     diffusion.simulate(dt);
 
     chem.getField().update(diffusion.getField());
 
-    // chem.getfield().update(diffusion.getfield());
+    // MSG("Chemistry start");
+    if (params.use_ai_surrogate) {
+      double ai_start_t = MPI_Wtime();
+      // Save current values from the tug field as predictor for the ai step
+      R["TMP"] = Rcpp::wrap(chem.getField().AsVector());
+      R.parseEval(
+          std::string("predictors <- setNames(data.frame(matrix(TMP, nrow=" +
+                      std::to_string(chem.getField().GetRequestedVecSize()) +
+                      ")), TMP_PROPS)"));
+      R.parseEval("predictors <- predictors[ai_surrogate_species]");
 
-    MSG("Chemistry step");
+      // Apply preprocessing
+      MSG("AI Preprocessing");
+      R.parseEval("predictors_scaled <- preprocess(predictors)");
+
+      // Predict
+      MSG("AI Prediction");
+      R.parseEval(
+          "aipreds_scaled <- prediction_step(model, predictors_scaled)");
+
+      // Apply postprocessing
+      MSG("AI Postprocessing");
+      R.parseEval("aipreds <- postprocess(aipreds_scaled)");
+
+      // Validate prediction and write valid predictions to chem field
+      MSG("AI Validation");
+      R.parseEval(
+          "validity_vector <- validate_predictions(predictors, aipreds)");
+
+      MSG("AI Marking accepted");
+      chem.set_ai_surrogate_validity_vector(R.parseEval("validity_vector"));
+
+      MSG("AI TempField");
+      std::vector<std::vector<double>> RTempField =
+          R.parseEval("set_valid_predictions(predictors,\
+                       aipreds,\
+                       validity_vector)");
+
+      MSG("AI Set Field");
+      Field predictions_field =
+          Field(R.parseEval("nrow(predictors)"), RTempField,
+                R.parseEval("colnames(predictors)"));
+
+      MSG("AI Update");
+      chem.getField().update(predictions_field);
+      double ai_end_t = MPI_Wtime();
+      R["ai_prediction_time"] = ai_end_t - ai_start_t;
+    }
 
     chem.simulate(dt);
 
+    /* AI surrogate iterative training*/
+    if (params.use_ai_surrogate) {
+      double ai_start_t = MPI_Wtime();
+
+      R["TMP"] = Rcpp::wrap(chem.getField().AsVector());
+      R.parseEval(
+          std::string("targets <- setNames(data.frame(matrix(TMP, nrow=" +
+                      std::to_string(chem.getField().GetRequestedVecSize()) +
+                      ")), TMP_PROPS)"));
+      R.parseEval("targets <- targets[ai_surrogate_species]");
+
+      // TODO: Check how to get the correct columns
+      R.parseEval("target_scaled <- preprocess(targets)");
+
+      MSG("AI: incremental training");
+      R.parseEval("model <- training_step(model, predictors_scaled, "
+                  "target_scaled, validity_vector)");
+      double ai_end_t = MPI_Wtime();
+      R["ai_training_time"] = ai_end_t - ai_start_t;
+    }
+
+    // MPI_Barrier(MPI_COMM_WORLD);
+    double end_t = MPI_Wtime();
+    dSimTime += end_t - start_t;
+    R["totaltime"] = dSimTime;
+
     // 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
@@ -289,13 +397,11 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, const RuntimeParameters &params,
 
     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;
+    // MSG();
   } // END SIMULATION LOOP
 
+  std::cout << std::endl;
+
   Rcpp::List chem_profiling;
   chem_profiling["simtime"] = chem.GetChemistryTime();
   chem_profiling["loop"] = chem.GetMasterLoopTime();
@@ -339,6 +445,97 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, const RuntimeParameters &params,
   return profiling;
 }
 
+std::vector<std::string> getSpeciesNames(const Field &&field, int root,
+                                         MPI_Comm comm) {
+  std::uint32_t n_elements;
+  std::uint32_t n_string_size;
+
+  int rank;
+  MPI_Comm_rank(comm, &rank);
+
+  const bool is_master = root == rank;
+
+  // first, the master sends all the species names iterative
+  if (is_master) {
+    n_elements = field.GetProps().size();
+    MPI_Bcast(&n_elements, 1, MPI_UINT32_T, root, MPI_COMM_WORLD);
+
+    for (std::uint32_t i = 0; i < n_elements; i++) {
+      n_string_size = field.GetProps()[i].size();
+      MPI_Bcast(&n_string_size, 1, MPI_UINT32_T, root, MPI_COMM_WORLD);
+      MPI_Bcast(const_cast<char *>(field.GetProps()[i].c_str()), n_string_size,
+                MPI_CHAR, root, MPI_COMM_WORLD);
+    }
+
+    return field.GetProps();
+  }
+
+  // now all the worker stuff
+  MPI_Bcast(&n_elements, 1, MPI_UINT32_T, root, comm);
+
+  std::vector<std::string> species_names_out(n_elements);
+
+  for (std::uint32_t i = 0; i < n_elements; i++) {
+    MPI_Bcast(&n_string_size, 1, MPI_UINT32_T, root, MPI_COMM_WORLD);
+
+    char recv_buf[n_string_size];
+
+    MPI_Bcast(recv_buf, n_string_size, MPI_CHAR, root, MPI_COMM_WORLD);
+
+    species_names_out[i] = std::string(recv_buf, n_string_size);
+  }
+
+  return species_names_out;
+}
+
+std::array<double, 2> getBaseTotals(Field &&field, int root, MPI_Comm comm) {
+  std::array<double, 2> base_totals;
+
+  int rank;
+  MPI_Comm_rank(comm, &rank);
+
+  const bool is_master = root == rank;
+
+  if (is_master) {
+    const auto h_col = field["H"];
+    const auto o_col = field["O"];
+
+    base_totals[0] = *std::min_element(h_col.begin(), h_col.end());
+    base_totals[1] = *std::min_element(o_col.begin(), o_col.end());
+    MPI_Bcast(base_totals.data(), 2, MPI_DOUBLE, root, MPI_COMM_WORLD);
+    return base_totals;
+  }
+
+  MPI_Bcast(base_totals.data(), 2, MPI_DOUBLE, root, comm);
+
+  return base_totals;
+}
+
+bool getHasID(Field &&field, int root, MPI_Comm comm) {
+  bool has_id;
+
+  int rank;
+  MPI_Comm_rank(comm, &rank);
+
+  const bool is_master = root == rank;
+
+  if (is_master) {
+    const auto ID_field = field["ID"];
+
+    std::set<double> unique_IDs(ID_field.begin(), ID_field.end());
+
+    has_id = unique_IDs.size() > 1;
+
+    MPI_Bcast(&has_id, 1, MPI_C_BOOL, root, MPI_COMM_WORLD);
+
+    return has_id;
+  }
+
+  MPI_Bcast(&has_id, 1, MPI_C_BOOL, root, comm);
+
+  return has_id;
+}
+
 int main(int argc, char *argv[]) {
   int world_size;
 
@@ -354,17 +551,24 @@ int main(int argc, char *argv[]) {
       MSG("Running POET version " + std::string(poet_version));
     }
 
+    init_global_functions(R);
+
     RuntimeParameters run_params;
 
-    switch (parseInitValues(argv, run_params)) {
-    case ParseRet::PARSER_ERROR:
-    case ParseRet::PARSER_HELP:
+    if (parseInitValues(argc, argv, run_params) != 0) {
       MPI_Finalize();
       return 0;
-    case ParseRet::PARSER_OK:
-      break;
     }
 
+    // switch (parseInitValues(argc, argv, run_params)) {
+    // case ParseRet::PARSER_ERROR:
+    // case ParseRet::PARSER_HELP:
+    //   MPI_Finalize();
+    //   return 0;
+    // case ParseRet::PARSER_OK:
+    //   break;
+    // }
+
     InitialList init_list(R);
     init_list.importList(run_params.init_params, MY_RANK != 0);
 
@@ -378,32 +582,54 @@ int main(int argc, char *argv[]) {
                               init_list.getChemistryInit(), MPI_COMM_WORLD);
 
     const ChemistryModule::SurrogateSetup surr_setup = {
-        init_list.getInitialGrid().GetProps(),
+        getSpeciesNames(init_list.getInitialGrid(), 0, MPI_COMM_WORLD),
+        getBaseTotals(init_list.getInitialGrid(), 0, MPI_COMM_WORLD),
+        getHasID(init_list.getInitialGrid(), 0, MPI_COMM_WORLD),
         run_params.use_dht,
         run_params.dht_size,
+        run_params.dht_snaps,
+        run_params.out_dir,
         run_params.use_interp,
         run_params.interp_bucket_entries,
         run_params.interp_size,
-        run_params.interp_min_entries};
+        run_params.interp_min_entries,
+        run_params.use_ai_surrogate};
 
     chemistry.masterEnableSurrogates(surr_setup);
 
     if (MY_RANK > 0) {
       chemistry.WorkerLoop();
     } else {
-
-      init_global_functions(R);
-
       // R.parseEvalQ("mysetup <- setup");
       // // if (MY_RANK == 0) { // get timestep vector from
       // // grid_init function ... //
-      *global_rt_setup =
-          master_init_R.value()(*global_rt_setup, run_params.out_dir,
-                                init_list.getInitialGrid().asSEXP());
+
+      *global_rt_setup = master_init_R(*global_rt_setup, run_params.out_dir,
+                                       init_list.getInitialGrid().asSEXP());
 
       // MDL: store all parameters
       // MSG("Calling R Function to store calling parameters");
       // R.parseEvalQ("StoreSetup(setup=mysetup)");
+      R["out_ext"] = run_params.out_ext;
+      R["out_dir"] = run_params.out_dir;
+
+      if (run_params.use_ai_surrogate) {
+        /* Incorporate ai surrogate from R */
+        R.parseEvalQ(ai_surrogate_r_library);
+        /* Use dht species for model input and output */
+        R["ai_surrogate_species"] =
+            init_list.getChemistryInit().dht_species.getNames();
+
+        const std::string ai_surrogate_input_script =
+            init_list.getChemistryInit().ai_surrogate_input_script;
+
+        MSG("AI: sourcing user-provided script");
+        R.parseEvalQ(ai_surrogate_input_script);
+
+        MSG("AI: initialize AI model");
+        R.parseEval("model <- initiate_model()");
+        R.parseEval("gpu_info()");
+      }
 
       MSG("Init done on process with rank " + std::to_string(MY_RANK));
 
@@ -420,14 +646,15 @@ int main(int argc, char *argv[]) {
 
       R["profiling"] = profiling;
       R["setup"] = *global_rt_setup;
+      R["setup$out_ext"] = run_params.out_ext;
 
-      string r_vis_code;
-      r_vis_code =
-          "saveRDS(profiling, file=paste0(setup$out_dir,'/timings.rds'));";
+      std::string r_vis_code;
+      r_vis_code = "SaveRObj(x = profiling, path = paste0(out_dir, "
+                   "'/timings.', setup$out_ext));";
       R.parseEval(r_vis_code);
 
       MSG("Done! Results are stored as R objects into <" + run_params.out_dir +
-          "/timings.rds>");
+          "/timings." + run_params.out_ext);
     }
   }
 

src/poet.hpp.in

@@ -23,7 +23,6 @@
 #pragma once
 
 #include <cstdint>
-#include <set>
 #include <string>
 #include <vector>
 
@@ -35,61 +34,39 @@ static const char *poet_version = "@POET_VERSION@";
 static const inline std::string kin_r_library = R"(@R_KIN_LIB@)";
 
 static const inline std::string init_r_library = R"(@R_INIT_LIB@)";
-
+static const inline std::string ai_surrogate_r_library =
+    R"(@R_AI_SURROGATE_LIB@)";
 static const inline std::string r_runtime_parameters = "mysetup";
 
-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"};
-
-constexpr uint32_t CHEM_DEFAULT_WORK_PACKAGE_SIZE = 32;
-
-constexpr uint32_t CHEM_DHT_SIZE_PER_PROCESS_MB = 1.5E3;
-
 struct RuntimeParameters {
   std::string out_dir;
   std::vector<double> timesteps;
 
-  bool print_progressbar;
-  uint32_t work_package_size;
-
   Rcpp::List init_params;
 
-  bool use_dht;
-  std::uint32_t dht_size;
-  std::uint8_t dht_snaps;
+  // MDL added to accomodate for qs::qsave/qread
+  bool as_rds = false;
+  bool as_qs = false;
+  std::string out_ext;
 
-  bool use_interp;
-  std::uint32_t interp_size;
-  std::uint32_t interp_min_entries;
-  std::uint32_t interp_bucket_entries;
+  bool print_progress = false;
 
-  struct ChemistryParams {
-    // std::string database_path;
-    // std::string input_script;
+  static constexpr std::uint32_t WORK_PACKAGE_SIZE_DEFAULT = 32;
+  std::uint32_t work_package_size = WORK_PACKAGE_SIZE_DEFAULT;
 
-    //   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_dht = false;
+  static constexpr std::uint32_t DHT_SIZE_DEFAULT = 1.5E3;
+  std::uint32_t dht_size = DHT_SIZE_DEFAULT;
+  static constexpr std::uint8_t DHT_SNAPS_DEFAULT = 0;
+  std::uint8_t dht_snaps = DHT_SNAPS_DEFAULT;
 
-    //   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;
+  bool use_interp = false;
+  static constexpr std::uint32_t INTERP_SIZE_DEFAULT = 100;
+  std::uint32_t interp_size = INTERP_SIZE_DEFAULT;
+  static constexpr std::uint32_t INTERP_MIN_ENTRIES_DEFAULT = 5;
+  std::uint32_t interp_min_entries = INTERP_MIN_ENTRIES_DEFAULT;
+  static constexpr std::uint32_t INTERP_BUCKET_ENTRIES_DEFAULT = 20;
+  std::uint32_t interp_bucket_entries = INTERP_BUCKET_ENTRIES_DEFAULT;
 
-    //   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);
-  };
+  bool use_ai_surrogate = false;
 };

test/testField.cpp

@@ -89,14 +89,14 @@ TEST_CASE("Field") {
   }
 
   SUBCASE("Apply R function (set Na to zero)") {
-    RHookFunction<Field> to_call(R, "simple_field");
+    poet::DEFunc to_call("simple_field");
     Field field_proc = to_call(dut.asSEXP());
 
     CHECK_EQ(field_proc["Na"], FieldColumn(dut.GetRequestedVecSize(), 0));
   }
 
   SUBCASE("Apply R function (add two fields)") {
-    RHookFunction<Field> to_call(R, "extended_field");
+    poet::DEFunc to_call("extended_field");
     Field field_proc = to_call(dut.asSEXP(), dut.asSEXP());
 
     CHECK_EQ(field_proc["Na"],

test/testNamedVector.cpp

@@ -9,7 +9,7 @@
 #include "testDataStructures.hpp"
 
 TEST_CASE("NamedVector") {
-  RInsidePOET &R = RInsidePOET::getInstance();
+  poet::RInsidePOET &R = poet::RInsidePOET::getInstance();
 
   R["sourcefile"] = RInside_source_file;
   R.parseEval("source(sourcefile)");
@@ -36,14 +36,14 @@ TEST_CASE("NamedVector") {
   }
 
   SUBCASE("Apply R function (set to zero)") {
-    RHookFunction<poet::NamedVector<double>> to_call(R, "simple_named_vec");
+    poet::DEFunc to_call("simple_named_vec");
     nv = to_call(nv);
 
     CHECK_EQ(nv[2], 0);
   }
 
   SUBCASE("Apply R function (second NamedVector)") {
-    RHookFunction<poet::NamedVector<double>> to_call(R, "extended_named_vec");
+    poet::DEFunc to_call("extended_named_vec");
 
     const std::vector<std::string> names{{"C", "H", "Mg"}};
     const std::vector<double> values{{0, 1, 2}};
@@ -56,8 +56,8 @@ TEST_CASE("NamedVector") {
   }
 
   SUBCASE("Apply R function (check if zero)") {
-    RHookFunction<bool> to_call(R, "bool_named_vec");
+    poet::DEFunc to_call("bool_named_vec");
 
-    CHECK_FALSE(to_call(nv));
+    CHECK_FALSE(Rcpp::as<bool>(to_call(nv)));
   }
 }