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 \
+ENV DEBIAN_FRONTEND=noninteractive
-    && sudo apt-get install -y \
+
-        cmake \
+RUN apt-get update \
+    && apt-get install -y \
+    sudo \
     git \
-        libeigen3-dev \
+    ninja-build \
-        libopenmpi-dev \
+    libmpfr-dev \
-        r-base-dev 
+    python3-dev && \
+    apt-get clean && \
+    rm -rf /var/lib/apt/lists/*
 
-RUN git clone https://github.com/doctest/doctest.git /doctest \
+WORKDIR /tmp
-    && cd /doctest \
-    && mkdir build \
-    && cd build \
-    && cmake .. \
-    && make install \
-    && cd / \
-    && rm -rf /doctest
 
-RUN  /usr/bin/R -q -e "install.packages(c('Rcpp', 'RInside'))"
+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.
+test:       # This job runs in the build stage, which runs first.
-  stage: build
+  stage: test
   script:
     - mkdir -p build && cd build
-    - cmake -DPOET_ENABLE_TESTING=ON -DPOET_PREPROCESS_BENCHS=OFF ..
+    - cmake -DPOET_ENABLE_TESTING=ON -DPOET_PREPROCESS_BENCHS=OFF -DCMAKE_BUILD_TYPE=Release ..
-    - make -j$(nproc)
-  cache: 
-    key: build-$CI_COMMIT_BRANCH
-    paths:
-      - build  
-  artifacts:
-    paths:
-      - build
-
-test-poet:
-  stage: test
-  dependencies:
-    - build-poet
-  script:
-    - cd build
     - 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"]
+[submodule "ext/litephreeqc"]
-	path = ext/iphreeqc
+	path = ext/litephreeqc
-	url = ../iphreeqc.git
+	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,8 +1,9 @@
 # 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
@@ -11,18 +12,16 @@ if(R_EXE)
 
 find_path(
   R_INCLUDE_DIR R.h
-    HINTS ${R_ROOT_DIR}
+  HINTS /usr/include /usr/local/include /usr/share ${R_ROOT_DIR}/include
-    PATHS /usr/include /usr/local/include /usr/share
+  PATH_SUFFIXES R/include R
-    PATH_SUFFIXES include/R R/include
+  REQUIRED
 )
 
 find_library(
   R_LIBRARY libR.so
   HINTS ${R_ROOT_DIR}/lib
+  REQUIRED
 )
-else()
-  message(FATAL_ERROR "No R runtime found!")
-endif()
 
 set(R_LIBRARIES ${R_LIBRARY})
 set(R_INCLUDE_DIRS ${R_INCLUDE_DIR})
@@ -71,11 +70,8 @@ list(APPEND R_INCLUDE_DIRS ${R_RInside_INCLUDE_DIR})
 # putting all together into interface library
 
 add_library(RRuntime INTERFACE IMPORTED)
-set_target_properties(
+target_link_libraries(RRuntime INTERFACE ${R_LIBRARIES})
-  RRuntime PROPERTIES
+target_include_directories(RRuntime INTERFACE ${R_INCLUDE_DIRS})
-  INTERFACE_LINK_LIBRARIES "${R_LIBRARIES}"
-  INTERFACE_INCLUDE_DIRECTORIES "${R_INCLUDE_DIRS}"
-)
 
 unset(R_LIBRARIES)
 unset(R_INCLUDE_DIRS)

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)
 

LICENSE

@@ -1,8 +1,8 @@
                     GNU GENERAL PUBLIC LICENSE
                        Version 2, June 1991
 
- Copyright (C) 1989, 1991 Free Software Foundation, Inc., <http://fsf.org/>
+ Copyright (C) 1989, 1991 Free Software Foundation, Inc.,
- 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ <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,25 +1,30 @@
 # POET
 
+**NOTE: GFZ is migrating its domain from <gfz-potsdam.de> to <gfz.de>.
+This should be finalized by the end of 2025. We adopt the NEW domain
+in all the links given below. If you encounter 'unreachable address'
+try the OLD domain.**
+
 [POET](https://doi.org/10.5281/zenodo.4757913) is a coupled reactive
 transport simulator implementing a parallel architecture and a fast,
 original MPI-based Distributed Hash Table.
 
-![POET's Coupling Scheme](./docs/Scheme_POET_en.svg)
+![POET's Coupling Scheme](./docs/POET_scheme.svg)
 
 ## Parsed code documentiation
 
 A parsed version of POET's documentation can be found at [Gitlab
-pages](https://naaice.git-pages.gfz-potsdam.de/poet).
+pages](https://naaice.git-pages.gfz.de/poet).
 
 ## External Libraries
 
 The following external libraries are shipped with POET:
 
 - **CLI11** - <https://github.com/CLIUtils/CLI11>
-- **IPhreeqc** with patches from GFZ/UP -
+- **litephreeqc**: IPhreeqc
-  <https://github.com/usgs-coupled/iphreeqc> -
+   (<https://github.com/usgs-coupled/iphreeqc>) with patches from
-  <https://git.gfz-potsdam.de/naaice/iphreeqc>
+   GFZ/UP: <https://git.gfz.de/naaice/litephreeqc>
-- **tug** - <https://git.gfz-potsdam.de/naaice/tug>
+- **tug** - <https://git.gfz.de/naaice/tug>
 
 ## Installation
 
@@ -29,7 +34,7 @@ To compile POET you need following software to be installed:
 
 - C/C++ compiler (tested with GCC)
 - MPI-Implementation (tested with OpenMPI and MVAPICH)
-- CMake 3.9+
+- CMake 3.20+
 - Eigen3 3.4+ (required by `tug`)
 - *optional*: `doxygen` with `dot` bindings for documentation
 - R language and environment including headers or `-dev` packages
@@ -41,7 +46,7 @@ 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
@@ -49,7 +54,7 @@ This can be simply achieved by issuing the following commands:
 $ R
 
 # install R dependencies (case sensitive!)
-> install.packages(c("Rcpp", "RInside","qs"))
+> install.packages(c("Rcpp", "RInside","qs","qs2"))
 > q(save="no")
 ```
 
@@ -59,7 +64,7 @@ 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-potsdam.de/naaice/poet.git
+git clone --recurse-submodules https://git.gfz.de/naaice/poet.git
 ```
 The `--recurse-submodules` option is a shorthand for:
 ```sh
@@ -74,7 +79,7 @@ usual:
 
 ```sh
 mkdir build && cd build
-cmake ..
+cmake -DCMAKE_BUILD_TYPE=Release ..
 ```
 
 This will create the directory `build` and processes the CMake files
@@ -110,7 +115,7 @@ follows:
 $ R
 
 # install R dependencies
-> install.packages(c("Rcpp", "RInside","qs"))
+> install.packages(c("Rcpp", "RInside","qs","qs2"))
 > q(save="no")
 
 # cd into POET project root
@@ -118,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
 ```
@@ -138,17 +143,17 @@ 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
 ```
 
@@ -182,7 +187,8 @@ The following parameters can be set:
 | **-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 RDS (.rds)    |
+| **--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                                       |
@@ -213,7 +219,7 @@ 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
@@ -226,7 +232,7 @@ 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
@@ -273,7 +279,7 @@ cp <project_root_dir>/bench/barite/{barite_50ai*,db_barite.dat,barite.pqi} .
 ./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.rds output
+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
@@ -284,7 +290,7 @@ produce any valid predictions.
 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).
+[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.
@@ -298,7 +304,7 @@ issue tracker or E-Mail.
 where: 
 
 - **output** - name of the output file (defaults to the input file
-  name with the extension `.rds`)
+  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

R_lib/init_r_lib.R

@@ -13,6 +13,34 @@
 ### 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
@@ -33,7 +61,7 @@ pqc_to_grid <- function(pqc_mat, grid) {
     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

R_lib/kin_r_library.R

@@ -94,7 +94,7 @@ master_iteration_end <- function(setup, state_T, state_C) {
     ## 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)
 }
@@ -109,69 +109,6 @@ 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)))
@@ -179,7 +116,7 @@ GetWorkPackageSizesVector <- function(n_packages, package_size, len) {
 
 
 ## Handler to read R objs from binary files using either builtin
-## readRDS() or qs::qread() based on file extension
+## readRDS(), qs::qread() or qs2::qs_read() based on file extension
 ReadRObj <- function(path) {
     ## code borrowed from tools::file_ext()
     pos <- regexpr("\\.([[:alnum:]]+)$", path)
@@ -187,20 +124,88 @@ ReadRObj <- function(path) {
 
     switch(extension,
         rds = readRDS(path),
-        qs  = qs::qread(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)
+    ## 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)
+        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

@@ -7,7 +7,7 @@ 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})
-        set(OUT_FILE_EXT ${OUT_FILE}.qs)
+        set(OUT_FILE_EXT ${OUT_FILE}.qs2)
 
         add_custom_command(
             OUTPUT ${OUT_FILE_EXT}

bench/barite/barite_200.R

@@ -35,15 +35,15 @@ diffusion_setup <- list(
 )
 
 dht_species <- c(
-  "H"         = 7,
+  "H"         = 3,
-  "O"         = 7,
+  "O"         = 3,
-  "Charge"    = 4,
+  "Charge"    = 6,
-  "Ba"        = 7,
+  "Ba"        = 6,
-  "Cl"        = 7,
+  "Cl"        = 6,
-  "S(6)"      = 7,
+  "S"         = 6,
-  "Sr"        = 7,
+  "Sr"        = 6,
-  "Barite"    = 4,
+  "Barite"    = 5,
-  "Celestite" = 4
+  "Celestite" = 5
 )
 
 chemistry_setup <- list(

bench/barite/barite_50ai.R

@@ -36,15 +36,15 @@ diffusion_setup <- list(
 )
 
 dht_species <- c(
-  "H"         = 4,
+  "H"         = 3,
-  "O"         = 10,
+  "O"         = 3,
-  "Charge"    = 4,
+  "Charge"    = 3,
-  "Ba"        = 7,
+  "Ba"        = 6,
-  "Cl"        = 4,
+  "Cl"        = 6,
-  "S(6)"      = 7,
+  "S"         = 6,
-  "Sr"        = 4,
+  "Sr"        = 6,
-  "Barite"    = 2,
+  "Barite"    = 5,
-  "Celestite" = 2
+  "Celestite" = 5
 )
 
 chemistry_setup <- list(

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)) {
+    # if ((old["Calcite"] == 0) != (new["Calcite"] == 0)) {
-        return(TRUE)
+    #     return(TRUE)
-    }
+    # }
-    if ((old["Dolomite"] == 0) != (new["Dolomite"] == 0)) {
+    # if ((old["Dolomite"] == 0) != (new["Dolomite"] == 0)) {
-        return(TRUE)
+    #     return(TRUE)
-    }
+    # }
     return(FALSE)
 }
 
-fuzz_input_dht_keys <- function(input) {
+# fuzz_input_dht_keys <- function(input) {
-    dht_species <- c(
+#     dht_species <- c(
-        "H" = 3,
+#         "H" = 3,
-        "O" = 3,
+#         "O" = 3,
-        "Charge" = 3,
+#         "Charge" = 3,
-        "C(4)" = 6,
+#         "C" = 6,
-        "Ca" = 6,
+#         "Ca" = 6,
-        "Cl" = 3,
+#         "Cl" = 3,
-        "Mg" = 5,
+#         "Mg" = 5,
-        "Calcite" = 4,
+#         "Calcite" = 4,
-        "Dolomite" = 4
+#         "Dolomite" = 4
-    )
+#     )
-    return(input[names(dht_species)])
+#     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/CMakeLists.txt

@@ -1,7 +1,4 @@
-add_library(POETLib)
+add_library(POETLib
-
-target_sources(POETLib 
-  PRIVATE
   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 
@@ -32,7 +38,7 @@ FetchContent_Declare(
   cli11
   QUIET 
   GIT_REPOSITORY https://github.com/CLIUtils/CLI11.git
-  GIT_TAG v2.4.2
+  GIT_TAG v2.5.0
 )
 
 FetchContent_MakeAvailable(cli11)

src/Chemistry/ChemistryDefs.hpp

@@ -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) {
+  WorkPackage(std::size_t _size)
-    input.resize(size);
+      : size(_size), input(size), output(size), mapping(size, CHEM_PQC) {}
-    output.resize(size);
-    mapping.resize(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) {
+    PhreeqcMatrix pqc_mat =
-      this->phreeqc_instances[pqc_id] =
+        PhreeqcMatrix(chem_params.database, chem_params.pqc_script,
-          std::make_unique<PhreeqcEngine>(pqc_config);
+                      chem_params.with_h0_o0, chem_params.with_redox);
-    }
+
-    // for (std::size_t i = 0; i < chem_params.pqc_ids.size(); i++) {
+    this->pqc_runner =
-    //   this->phreeqc_instances[chem_params.pqc_ids[i]] =
+        std::make_unique<PhreeqcRunner>(pqc_mat.subset(chem_params.pqc_ids));
-    //       std::make_unique<PhreeqcWrapper>(
-    //           chem_params.database, chem_params.pqc_scripts[i],
-    //           chem_params.pqc_sol_order, chem_params.field_header, wp_size_);
-    // }
   }
 }
 
@@ -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,7 +257,8 @@ 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
+            static_cast<std::uint32_t>(map_copy[i]) -
+            (map_copy[i] > InterpolationModule::COARSE_DIFF
                  ? InterpolationModule::COARSE_DIFF
                  : 0);
         map_copy[i] = signif;
@@ -372,6 +318,7 @@ void poet::ChemistryModule::unshuffleField(const std::vector<double> &in_buffer,
   }
 }
 
-void poet::ChemistryModule::set_ai_surrogate_validity_vector(std::vector<int> r_vector) {
+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,9 +76,13 @@ 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;
@@ -95,8 +100,15 @@ public:
     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) {
@@ -239,7 +251,8 @@ public:
 
 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);
 
@@ -390,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/SurrogateModels/DHT.c

@@ -105,6 +105,8 @@ DHT *DHT_create(MPI_Comm comm, uint64_t size, unsigned int data_size,
   object->index_count = 9 - (index_bytes / 8);
   object->index = (uint64_t *)malloc((object->index_count) * sizeof(uint64_t));
   object->mem_alloc = mem_alloc;
+  object->sum_idx = 0;
+  object->cnt_idx = 0;
 
   // if set, initialize dht_stats
 #ifdef DHT_STATISTICS
@@ -189,6 +191,9 @@ int DHT_write_accumulate(DHT *table, const void *send_key, int data_size,
     }
   }
 
+  table->cnt_idx += 1;
+  table->sum_idx += (i + 1);
+
   if (result == DHT_WRITE_SUCCESS_WITH_EVICTION) {
     memset((char *)table->send_entry + 1 + table->key_size, '\0',
            table->data_size);
@@ -279,6 +284,9 @@ int DHT_write(DHT *table, void *send_key, void *send_data, uint32_t *proc,
     }
   }
 
+  table->cnt_idx += 1;
+  table->sum_idx += (i + 1);
+
   // put data to DHT (with last selected index by value i)
   if (MPI_Put(table->send_entry, 1 + table->data_size + table->key_size,
               MPI_BYTE, dest_rank, table->index[i],
@@ -550,6 +558,41 @@ int DHT_free(DHT *table, int *eviction_counter, int *readerror_counter) {
   return DHT_SUCCESS;
 }
 
+float DHT_get_used_idx_factor(DHT *table, int with_reset) {
+  int rank;
+  MPI_Comm_rank(table->communicator, &rank);
+
+  float my_avg_idx = (float)table->sum_idx / (float)table->cnt_idx;
+
+  float max_mean_index;
+
+  MPI_Reduce(&my_avg_idx, &max_mean_index, 1, MPI_FLOAT, MPI_MAX, 0,
+             table->communicator);
+
+  MPI_Bcast(&max_mean_index, 1, MPI_FLOAT, 0, table->communicator);
+
+  if (!!with_reset) {
+    table->sum_idx = 0;
+    table->cnt_idx = 0;
+  }
+
+  return max_mean_index;
+}
+
+int DHT_flush(DHT *table) {
+  // make sure all processes are synchronized
+  MPI_Barrier(table->communicator);
+
+  // wipe local memory with zeros
+  memset(table->mem_alloc, '\0',
+         table->table_size * (1 + table->data_size + table->key_size));
+
+  table->sum_idx = 0;
+  table->cnt_idx = 0;
+
+  return DHT_SUCCESS;
+}
+
 int DHT_print_statistics(DHT *table) {
 #ifdef DHT_STATISTICS
   int *written_buckets;

src/Chemistry/SurrogateModels/DHT.h

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

src/Chemistry/SurrogateModels/DHT_Wrapper.cpp

@@ -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,7 +129,13 @@ 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) {
+    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!
@@ -137,7 +144,7 @@ void DHT_Wrapper::fillDHT(const WorkPackage &work_package) {
       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)) {
+      if (Rcpp::as<bool>(hooks.dht_fill(old_values, new_values))) {
         continue;
       }
     }
@@ -165,7 +172,6 @@ void DHT_Wrapper::fillDHT(const WorkPackage &work_package) {
     dht_results.filledDHT[i] = true;
   }
 }
-}
 
 inline std::vector<double>
 DHT_Wrapper::outputToInputAndRates(const std::vector<double> &old_results,
@@ -270,7 +276,7 @@ LookupKey DHT_Wrapper::fuzzForDHT_R(const std::vector<double> &cell,
   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;
 
@@ -290,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;
 }
@@ -297,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;
@@ -323,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

@@ -87,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
    *
@@ -264,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/Interpolation.hpp

@@ -166,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;
@@ -261,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

@@ -14,6 +14,7 @@
 #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;
@@ -116,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

@@ -10,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 {
@@ -23,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/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>
@@ -48,13 +47,15 @@ void poet::ChemistryModule::WorkerLoop() {
     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?
+        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);
+      MPI_Bcast(&this->ai_surrogate_validity_vector.front(), this->n_cells,
+                MPI_INT, 0, this->group_comm);
       break;
     }
     case CHEM_WORK_LOOP: {
@@ -187,7 +188,6 @@ void poet::ChemistryModule::WorkerDoWork(MPI_Status &probe_status,
     }
   }
 
-
   phreeqc_time_start = MPI_Wtime();
 
   WorkerRunWorkPackage(s_curr_wp, current_sim_time, dt);
@@ -246,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()");
@@ -300,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);
     }
-
-    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++];
   }
+  this->pqc_runner->run(inout_chem, dTimestep, to_ignore);
+
+  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"]);
   }
@@ -69,21 +67,16 @@ 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++) {
+  // for (std::size_t i = 0; i < pqc_ids.size(); ++i) {
-    POETInitCell cell = {
+  //   chem_init.pqc_input[pqc_ids[i]] = pqc_scripts[i];
-        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};
-  }
 
   // chem_init.pqc_sol_order = pqc_solutions;
 

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>
@@ -53,98 +52,102 @@ static std::vector<TugType> colMajToRowMaj(const Rcpp::NumericVector &vec,
   }
 }
 
-static std::vector<std::string>
+// static std::vector<std::string>
-extend_transport_names(std::unique_ptr<PhreeqcInit> &phreeqc,
+// extend_transport_names(std::unique_ptr<PhreeqcInit> &phreeqc,
-                       const Rcpp::List &boundaries_list,
+//                        const Rcpp::List &boundaries_list,
-                       const Rcpp::List &inner_boundaries,
+//                        const Rcpp::List &inner_boundaries,
-                       const std::vector<std::string> &old_trans_names) {
+//                        const std::vector<std::string> &old_trans_names) {
 
-  std::vector<std::string> transport_names = old_trans_names;
+//   std::vector<std::string> transport_names = old_trans_names;
-  std::set<int> constant_pqc_ids;
+//   std::set<int> constant_pqc_ids;
 
-  for (const auto &side : tug_side_mapping) {
+//   for (const auto &side : tug_side_mapping) {
-    if (!boundaries_list.containsElementNamed(side.second.c_str())) {
+//     if (!boundaries_list.containsElementNamed(side.second.c_str())) {
-      continue;
+//       continue;
-    }
+//     }
 
-    Rcpp::List mapping = boundaries_list[side.second];
+//     Rcpp::List mapping = boundaries_list[side.second];
 
-    const Rcpp::NumericVector cells = mapping["cell"];
+//     const Rcpp::NumericVector cells = mapping["cell"];
-    const Rcpp::NumericVector values = mapping["sol_id"];
+//     const Rcpp::NumericVector values = mapping["sol_id"];
-    const Rcpp::CharacterVector type_str = mapping["type"];
+//     const Rcpp::CharacterVector type_str = mapping["type"];
 
-    if (cells.size() != values.size()) {
+//     if (cells.size() != values.size()) {
-      throw std::runtime_error("Boundary [" + side.second +
+//       throw std::runtime_error("Boundary [" + side.second +
-                               "] cells and values are not the same "
+//                                "] cells and values are not the same "
-                               "length");
+//                                "length");
-    }
+//     }
 
-    for (auto i = 0; i < cells.size(); i++) {
+//     for (auto i = 0; i < cells.size(); i++) {
-      if (type_str[i] == "constant") {
+//       if (type_str[i] == "constant") {
-        constant_pqc_ids.insert(static_cast<std::size_t>(values[i]));
+//         constant_pqc_ids.insert(static_cast<std::size_t>(values[i]));
-      }
+//       }
-    }
+//     }
-  }
+//   }
 
-  if (inner_boundaries.size() > 0) {
+//   if (inner_boundaries.size() > 0) {
-    const Rcpp::NumericVector values = inner_boundaries["sol_id"];
+//     const Rcpp::NumericVector values = inner_boundaries["sol_id"];
-    for (auto i = 0; i < values.size(); i++) {
+//     for (auto i = 0; i < values.size(); i++) {
-      constant_pqc_ids.insert(static_cast<std::size_t>(values[i]));
+//       constant_pqc_ids.insert(static_cast<std::size_t>(values[i]));
-    }
+//     }
-  }
+//   }
 
-  if (!constant_pqc_ids.empty()) {
+//   if (!constant_pqc_ids.empty()) {
-    constexpr std::size_t keep_h_o_charge = 3;
+//     constexpr std::size_t keep_h_o_charge = 3;
 
-    for (const auto &pqc_id : constant_pqc_ids) {
+//     for (const auto &pqc_id : constant_pqc_ids) {
-      const auto solution_names = phreeqc->getSolutionNames(pqc_id);
+//       const auto solution_names = phreeqc->getSolutionNames(pqc_id);
 
-      // add those strings which are not already in the transport_names
+//       // add those strings which are not already in the transport_names
-      for (const auto &name : solution_names) {
+//       for (const auto &name : solution_names) {
-        if (std::find(transport_names.begin(), transport_names.end(), name) ==
+//         if (std::find(transport_names.begin(), transport_names.end(), name)
-            transport_names.end()) {
+//         ==
-          auto position =
+//             transport_names.end()) {
-              std::lower_bound(transport_names.begin() + keep_h_o_charge,
+//           auto position =
-                               transport_names.end(), name);
+//               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
+// static Rcpp::List
-extend_initial_grid(const Rcpp::List &initial_grid,
+// extend_initial_grid(const Rcpp::List &initial_grid,
-                    const std::vector<std::string> &transport_names) {
+//                     const std::vector<std::string> &transport_names) {
-  // std::vector<std::string> names_to_add(transport_names.begin() + old_size,
+//   // std::vector<std::string> names_to_add(transport_names.begin() +
-  //                                       transport_names.end());
+//   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->transport_names = extend_transport_names(
-      this->phreeqc, boundaries_list, inner_boundaries, this->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) {
+  // if (old_transport_size != new_transport_size) {
-    this->initial_grid =
+  //   this->initial_grid =
-        extend_initial_grid(this->initial_grid, this->transport_names);
+  //       extend_initial_grid(this->initial_grid, this->transport_names);
-  }
+  // }
 
   for (const auto &species_name : this->transport_names) {
     Rcpp::List spec_list;
@@ -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>(
+            c_value[c_id] =
-                resolve_R(this->phreeqc_mat, Rcpp::wrap(species_name), values[i]));
+                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");
           }
@@ -211,8 +217,11 @@ 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(
+        c_value.push_back(static_cast<TugType>(
-            this->phreeqc_mat, Rcpp::wrap(species_name), inner_pqc_id_vec[i])));
+            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_name] =
@@ -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;
 

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) {
+static Rcpp::List expandGrid(const PhreeqcMatrix &pqc_mat,
-  // parse input values
+                             const std::vector<int> unique_ids,
-  Rcpp::Function unique_R("unique");
+                             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);
+  bool with_h0_o0 =
-  this->initial_grid = matToGrid(R, this->phreeqc_mat, grid_def);
+      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);
+  if (with_h0_o0 && !with_redox) {
-
+    throw std::runtime_error(
-  std::vector<std::string> colnames =
+        "Output of H(0) and O(0) can only be used with redox.");
-      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));
   }
+
+  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

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]);
+  PhreeqcMatrix phreeqc = prepareGrid(setup[grid_key]);
-  initDiffusion(setup[diffusion_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>>(
+  //   this->pqc_solutions = Rcpp::as<std::vector<std::string>>(
-      setup[static_cast<int>(ExportList::CHEM_PQC_SOLUTIONS)]);
+  //       setup[static_cast<int>(ExportList::CHEM_PQC_SOLUTIONS)]);
-  this->pqc_solution_primaries = Rcpp::as<std::vector<std::string>>(
+  //   this->pqc_solution_primaries = Rcpp::as<std::vector<std::string>>(
-      setup[static_cast<int>(ExportList::CHEM_PQC_SOLUTION_PRIMARY)]);
+  //       setup[static_cast<int>(ExportList::CHEM_PQC_SOLUTION_PRIMARY)]);
-  this->pqc_exchanger =
+  //   this->pqc_exchanger =
-      Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_EXCHANGER)]);
+  //       Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_EXCHANGER)]);
-  this->pqc_kinetics =
+  //   this->pqc_kinetics =
-      Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_KINETICS)]);
+  //       Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_KINETICS)]);
-  this->pqc_equilibrium =
+  //   this->pqc_equilibrium =
-      Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_EQUILIBRIUM)]);
+  //       Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_EQUILIBRIUM)]);
-  this->pqc_surface_comps =
+  //   this->pqc_surface_comps =
-      Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_SURFACE_COMPS)]);
+  //       Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_SURFACE_COMPS)]);
-  this->pqc_surface_charges =
+  //   this->pqc_surface_charges =
-      Rcpp::List(setup[static_cast<int>(ExportList::CHEM_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)]);
@@ -108,25 +113,29 @@ 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)] =
+  //   out[static_cast<int>(ExportList::CHEM_PQC_SOLUTIONS)] =
-      Rcpp::wrap(this->pqc_solutions);
+  //       Rcpp::wrap(this->pqc_solutions);
-  out[static_cast<int>(ExportList::CHEM_PQC_SOLUTION_PRIMARY)] =
+  //   out[static_cast<int>(ExportList::CHEM_PQC_SOLUTION_PRIMARY)] =
-      Rcpp::wrap(this->pqc_solution_primaries);
+  //       Rcpp::wrap(this->pqc_solution_primaries);
-  out[static_cast<int>(ExportList::CHEM_PQC_EXCHANGER)] =
+  //   out[static_cast<int>(ExportList::CHEM_PQC_EXCHANGER)] =
-      Rcpp::wrap(this->pqc_exchanger);
+  //       Rcpp::wrap(this->pqc_exchanger);
-  out[static_cast<int>(ExportList::CHEM_PQC_KINETICS)] =
+  //   out[static_cast<int>(ExportList::CHEM_PQC_KINETICS)] =
-      Rcpp::wrap(this->pqc_kinetics);
+  //       Rcpp::wrap(this->pqc_kinetics);
-  out[static_cast<int>(ExportList::CHEM_PQC_EQUILIBRIUM)] =
+  //   out[static_cast<int>(ExportList::CHEM_PQC_EQUILIBRIUM)] =
-      Rcpp::wrap(this->pqc_equilibrium);
+  //       Rcpp::wrap(this->pqc_equilibrium);
-  out[static_cast<int>(ExportList::CHEM_PQC_SURFACE_COMPS)] =
+  //   out[static_cast<int>(ExportList::CHEM_PQC_SURFACE_COMPS)] =
-      Rcpp::wrap(this->pqc_surface_comps);
+  //       Rcpp::wrap(this->pqc_surface_comps);
-  out[static_cast<int>(ExportList::CHEM_PQC_SURFACE_CHARGES)] =
+  //   out[static_cast<int>(ExportList::CHEM_PQC_SURFACE_CHARGES)] =
-      Rcpp::wrap(this->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);

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,16 +51,11 @@ private:
     DIFFU_ALPHA_X,
     DIFFU_ALPHA_Y,
     CHEM_DATABASE,
-    CHEM_FIELD_HEADER,
+    CHEM_PQC_SCRIPT,
-    CHEM_PQC_SCRIPTS,
+    CHEM_PQC_WITH_H0_O0,
+    CHEM_PQC_WITH_REDOX,
     CHEM_PQC_IDS,
-    CHEM_PQC_SOLUTIONS,
+    CHEM_FIELD_HEADER,
-    CHEM_PQC_SOLUTION_PRIMARY,
-    CHEM_PQC_EXCHANGER,
-    CHEM_PQC_KINETICS,
-    CHEM_PQC_EQUILIBRIUM,
-    CHEM_PQC_SURFACE_COMPS,
-    CHEM_PQC_SURFACE_CHARGES,
     CHEM_DHT_SPECIES,
     CHEM_INTERP_SPECIES,
     CHEM_HOOKS,
@@ -76,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,
@@ -89,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",
+      GridMembersString = {"pqc_in_string", "pqc_in_file",   "pqc_with_redox",
-                           "pqc_db_file",    "grid_def",    "grid_size",
+                           "pqc_wth_h0_o0", "pqc_db_string", "pqc_db_file",
-                           "constant_cells", "porosity"};
+                           "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};
 
@@ -114,9 +110,6 @@ private:
 
   Rcpp::List initial_grid;
 
-  // No export
-  Rcpp::NumericMatrix phreeqc_mat;
-
 public:
   struct DiffusionInit {
     using BoundaryMap =
@@ -169,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;
@@ -190,17 +185,12 @@ 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;
@@ -227,10 +217,12 @@ public:
     // std::vector<std::string> field_header;
 
     std::string database;
-    // std::vector<std::string> pqc_scripts;
+    std::string pqc_script;
-    // std::vector<int> pqc_ids;
+    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;
 

src/Transport/DiffusionModule.cpp

@@ -1,24 +1,7 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
 /*
-** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
+ * DiffusionModule.cpp - Coupling module betweeen POET and tug
-** Potsdam)
+ * Copyright (C) 2018-2025 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.
  */
 
 #include "DiffusionModule.hpp"
@@ -26,72 +9,40 @@
 #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);
+    tug::RowMajMatMap<TugType> alpha_x_map(
-    Eigen::MatrixX<TugType> alpha_x =
+        this->param_list.alpha_x[sol_name].data(), n_rows, n_cols);
-        VecToMatrix(this->param_list.alpha_x[sol_name], n_rows, n_cols);
+    tug::RowMajMatMap<TugType> alpha_y_map(
-    Eigen::MatrixX<TugType> alpha_y =
+        this->param_list.alpha_y[sol_name].data(), n_rows, n_cols);
-        VecToMatrix(this->param_list.alpha_y[sol_name], 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);
+    diffusion_solver.setTimestep(requested_dt);
-    grid.setConcentrations(conc);
+    diffusion_solver.run();
-
-    sim.setTimestep(requested_dt);
-
-    sim.run();
-
-    species_conc = MatrixToVec(grid.getConcentrations());
   }
 
   const auto end_diffusion_t = std::chrono::high_resolution_clock::now();

src/initializer.cpp

@@ -34,9 +34,11 @@ int main(int argc, char **argv) {
                "input script")
       ->default_val(false);
 
-  bool asRDS;
+  bool asRDS{false};
-  app.add_flag("-r, --rds", asRDS, "Save output as .rds file instead of .qs")
+  app.add_flag("-r, --rds", asRDS, "Save output as .rds")->default_val(false);
-      ->default_val(false);
+
+  bool asQS{false};
+  app.add_flag("-q, --qs", asQS, "Save output as .qs")->default_val(false);
 
   CLI11_PARSE(app, argc, argv);
 
@@ -69,7 +71,13 @@ int main(int argc, char **argv) {
   }
 
   // append the correct file extension
-  output_file += asRDS ? ".rds" : ".qs";
+  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) {

src/poet.cpp

@@ -34,10 +34,13 @@
 #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 <set>
 #include <string>
 
 #include <CLI/CLI.hpp>
@@ -45,7 +48,7 @@
 #include <poet.hpp>
 #include <vector>
 
-using namespace std;
+// using namespace std;
 using namespace poet;
 using namespace Rcpp;
 
@@ -57,7 +60,7 @@ static std::unique_ptr<Rcpp::List> global_rt_setup;
 // before the R runtime is initialized
 static poet::DEFunc master_init_R;
 static poet::DEFunc master_iteration_end_R;
-static poet::DEFunc store_setup_R;
+// MDL: unused -> static poet::DEFunc store_setup_R;
 static poet::DEFunc ReadRObj_R;
 static poet::DEFunc SaveRObj_R;
 static poet::DEFunc source_R;
@@ -66,7 +69,7 @@ static void init_global_functions(RInside &R) {
   R.parseEval(kin_r_library);
   master_init_R = DEFunc("master_init");
   master_iteration_end_R = DEFunc("master_iteration_end");
-  store_setup_R = DEFunc("StoreSetup");
+  // MDL: unused -> store_setup_R = DEFunc("StoreSetup");
   source_R = DEFunc("source");
   ReadRObj_R = DEFunc("ReadRObj");
   SaveRObj_R = DEFunc("SaveRObj");
@@ -146,7 +149,10 @@ int parseInitValues(int argc, char **argv, RuntimeParameters &params) {
                "Enable AI surrogate for chemistry module");
 
   app.add_flag("--rds", params.as_rds,
-               "Save output as .rds file instead of .qs");
+               "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;
@@ -174,7 +180,11 @@ int parseInitValues(int argc, char **argv, RuntimeParameters &params) {
   }
 
   // set the output extension
-  params.out_ext = params.as_rds ? "rds" : "qs";
+  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));
@@ -287,18 +297,20 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, const RuntimeParameters &params,
 
     const double &dt = params.timesteps[iter - 1];
 
-    //  cout << "CPP: Next time step is " << dt << "[s]" << endl;
+    std::cout << std::endl;
-    MSG("Next time step is " + std::to_string(dt) + " [s]");
 
     /* displaying iteration number, with C++ and R iterator */
-    MSG("Going through iteration " + std::to_string(iter));
+    MSG("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());
 
-    MSG("Chemistry step");
+    // 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
@@ -314,16 +326,16 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, const RuntimeParameters &params,
       R.parseEval("predictors_scaled <- preprocess(predictors)");
 
       // Predict
-      MSG("AI Predict");
+      MSG("AI Prediction");
       R.parseEval(
           "aipreds_scaled <- prediction_step(model, predictors_scaled)");
 
       // Apply postprocessing
-      MSG("AI Postprocesing");
+      MSG("AI Postprocessing");
       R.parseEval("aipreds <- postprocess(aipreds_scaled)");
 
       // Validate prediction and write valid predictions to chem field
-      MSG("AI Validate");
+      MSG("AI Validation");
       R.parseEval(
           "validity_vector <- validate_predictions(predictors, aipreds)");
 
@@ -385,9 +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();
+    // MSG();
   } // END SIMULATION LOOP
 
+  std::cout << std::endl;
+
   Rcpp::List chem_profiling;
   chem_profiling["simtime"] = chem.GetChemistryTime();
   chem_profiling["loop"] = chem.GetMasterLoopTime();
@@ -474,6 +488,54 @@ std::vector<std::string> getSpeciesNames(const Field &&field, int root,
   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;
 
@@ -520,10 +582,13 @@ int main(int argc, char *argv[]) {
                               init_list.getChemistryInit(), MPI_COMM_WORLD);
 
     const ChemistryModule::SurrogateSetup surr_setup = {
-
         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,
@@ -583,7 +648,7 @@ int main(int argc, char *argv[]) {
       R["setup"] = *global_rt_setup;
       R["setup$out_ext"] = run_params.out_ext;
 
-      string r_vis_code;
+      std::string r_vis_code;
       r_vis_code = "SaveRObj(x = profiling, path = paste0(out_dir, "
                    "'/timings.', setup$out_ext));";
       R.parseEval(r_vis_code);

src/poet.hpp.in

@@ -23,7 +23,6 @@
 #pragma once
 
 #include <cstdint>
-#include <set>
 #include <string>
 #include <vector>
 
@@ -45,27 +44,29 @@ struct RuntimeParameters {
 
   Rcpp::List init_params;
 
+  // MDL added to accomodate for qs::qsave/qread
   bool as_rds = false;
-  std::string out_ext; // MDL added to accomodate for qs::qsave/qread
+  bool as_qs = false;
+  std::string out_ext;
 
   bool print_progress = false;
 
   static constexpr std::uint32_t WORK_PACKAGE_SIZE_DEFAULT = 32;
-  std::uint32_t work_package_size;
+  std::uint32_t work_package_size = WORK_PACKAGE_SIZE_DEFAULT;
 
   bool use_dht = false;
   static constexpr std::uint32_t DHT_SIZE_DEFAULT = 1.5E3;
-  std::uint32_t dht_size;
+  std::uint32_t dht_size = DHT_SIZE_DEFAULT;
   static constexpr std::uint8_t DHT_SNAPS_DEFAULT = 0;
-  std::uint8_t dht_snaps;
+  std::uint8_t dht_snaps = DHT_SNAPS_DEFAULT;
 
   bool use_interp = false;
   static constexpr std::uint32_t INTERP_SIZE_DEFAULT = 100;
-  std::uint32_t interp_size;
+  std::uint32_t interp_size = INTERP_SIZE_DEFAULT;
   static constexpr std::uint32_t INTERP_MIN_ENTRIES_DEFAULT = 5;
-  std::uint32_t interp_min_entries;
+  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;
+  std::uint32_t interp_bucket_entries = INTERP_BUCKET_ENTRIES_DEFAULT;
 
   bool use_ai_surrogate = false;
 };