Merge branch 'mdl/defaultstorelease' into 'main'

DEFAULT_BUILD_TYPE set to Release See merge request naaice/poet!60
set DEFAULT_BUILD_TYPE to "Release"; defined properties so one cycles through when using ccmake
2025-12-15 12:28:22 +01:00 · 2025-12-06 22:37:17 +01:00 · 2025-11-14 19:19:07 +01:00 · 2025-11-14 15:31:46 +01:00 · 2025-11-14 15:28:25 +01:00 · 2025-11-13 16:45:53 +01:00
49 changed files with 1727 additions and 764 deletions
--- a/.devcontainer/Dockerfile
+++ b/.devcontainer/Dockerfile
@ -1,20 +1,108 @@
-FROM mcr.microsoft.com/vscode/devcontainers/base:debian
+FROM gcc:11.2.0 AS builder

-RUN sudo apt-get update && export DEBIAN_FRONTEND=noninteractive \
-    && sudo apt-get install -y \
-        cmake \
-        git \
-        libeigen3-dev \
-        libopenmpi-dev \
-        r-base-dev 
+ENV DEBIAN_FRONTEND=noninteractive

-RUN git clone https://github.com/doctest/doctest.git /doctest \
-    && cd /doctest \
-    && mkdir build \
-    && cd build \
-    && cmake .. \
-    && make install \
-    && cd / \
-    && rm -rf /doctest
+RUN apt-get update \
+    && apt-get install -y \
+    sudo \
+    git \
+    ninja-build \
+    libmpfr-dev \
+    python3-dev && \
+    apt-get clean && \
+    rm -rf /var/lib/apt/lists/*

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

 stages:          # List of stages for jobs, and their order of execution
  - release
-  - build
  - test

 variables:
@ -28,26 +27,11 @@ variables:
  SOURCE_ARCHIVE_NAME: 'poet_${CI_COMMIT_TAG}_sources.tar.gz'
  CHANGELOG_FILE: 'commit_changelog.md'

-build-poet:       # This job runs in the build stage, which runs first.
-  stage: build
+test:       # This job runs in the build stage, which runs first.
+  stage: test
  script:
    - mkdir -p build && cd build
-    - cmake -DPOET_ENABLE_TESTING=ON -DPOET_PREPROCESS_BENCHS=OFF ..
-    - make -j$(nproc)
-  cache: 
-    key: build-$CI_COMMIT_BRANCH
-    paths:
-      - build  
-  artifacts:
-    paths:
-      - build
-
-test-poet:
-  stage: test
-  dependencies:
-    - build-poet
-  script:
-    - cd build
+    - cmake -DPOET_ENABLE_TESTING=ON -DPOET_PREPROCESS_BENCHS=OFF -DCMAKE_BUILD_TYPE=Release ..
    - make -j$(nproc) check

 pages:
@ -65,6 +49,22 @@ pages:
  rules:
    - if: $CI_COMMIT_REF_NAME == $CI_DEFAULT_BRANCH || $CI_COMMIT_TAG

+push:
+  stage: release
+  variables:
+    GITHUB_REPOSITORY: 'git@github.com:POET-Simulator/POET.git'
+  before_script:
+    # I know that there is this file env variable in gitlab, but somehow it does not work for me (still complaining about white spaces ...)
+    # Therefore, the ssh key is stored as a base64 encoded string
+    - mkdir -p ~/.ssh && echo $GITHUB_SSH_PRIVATE_KEY | base64 -d > ~/.ssh/id_ed25519 && chmod 0600 ~/.ssh/id_ed25519
+    - ssh-keyscan github.com >> ~/.ssh/known_hosts
+    - echo $MIRROR_SCRIPT | base64 -d > mirror.sh && chmod +x mirror.sh
+  script:
+    - if [[-d poet.git ]]; then rm -rf poet.git; fi
+    - git clone --mirror "https://git.gfz-potsdam.de/naaice/poet.git" "poet.git" && cd poet.git
+    - git push --mirror $GITHUB_REPOSITORY
+  allow_failure: true
+
 #archive-sources:       # This job runs in the build stage, which runs first.
 #  image: python:3
 #  stage: release
--- a/.gitmodules
+++ b/.gitmodules
@ -2,6 +2,6 @@
 	path = ext/tug
 	url = ../tug.git

-[submodule "ext/iphreeqc"]
-	path = ext/iphreeqc
-	url = ../iphreeqc.git
+[submodule "ext/litephreeqc"]
+	path = ext/litephreeqc
+	url = ../litephreeqc.git
--- a/CITATION.cff
+++ b/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
--- a/CMake/FindRRuntime.cmake
+++ b/CMake/FindRRuntime.cmake
@ -1,28 +1,27 @@
 # prepare R environment (Rcpp + RInside)
-find_program(R_EXE "R")
+find_program(R_EXE "R"
+  REQUIRED
+)

 # search for R executable, R header file and library path
-if(R_EXE)
-  execute_process(
-    COMMAND ${R_EXE} RHOME
-    OUTPUT_VARIABLE R_ROOT_DIR
-    OUTPUT_STRIP_TRAILING_WHITESPACE
-  )
+execute_process(
+  COMMAND ${R_EXE} RHOME
+  OUTPUT_VARIABLE R_ROOT_DIR
+  OUTPUT_STRIP_TRAILING_WHITESPACE
+)

-  find_path(
-    R_INCLUDE_DIR R.h
-    HINTS ${R_ROOT_DIR}
-    PATHS /usr/include /usr/local/include /usr/share
-    PATH_SUFFIXES include/R R/include
-  )
+find_path(
+  R_INCLUDE_DIR R.h
+  HINTS /usr/include /usr/local/include /usr/share ${R_ROOT_DIR}/include
+  PATH_SUFFIXES R/include R
+  REQUIRED
+)

-  find_library(
-    R_LIBRARY libR.so
-    HINTS ${R_ROOT_DIR}/lib
-  )
-else()
-  message(FATAL_ERROR "No R runtime found!")
-endif()
+find_library(
+  R_LIBRARY libR.so
+  HINTS ${R_ROOT_DIR}/lib
+  REQUIRED
+)

 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(
-  RRuntime PROPERTIES
-  INTERFACE_LINK_LIBRARIES "${R_LIBRARIES}"
-  INTERFACE_INCLUDE_DIRECTORIES "${R_INCLUDE_DIRS}"
-)
+target_link_libraries(RRuntime INTERFACE ${R_LIBRARIES})
+target_include_directories(RRuntime INTERFACE ${R_INCLUDE_DIRS})

 unset(R_LIBRARIES)
 unset(R_INCLUDE_DIRS)
--- a/CMakeLists.txt
+++ b/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)

--- a/62
+++ b/62
@ -1,8 +1,8 @@
                    GNU GENERAL PUBLIC LICENSE
                       Version 2, June 1991

- Copyright (C) 1989, 1991 Free Software Foundation, Inc., <http://fsf.org/>
- 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ Copyright (C) 1989, 1991 Free Software Foundation, Inc.,
+ <https://fsf.org/>
 Everyone is permitted to copy and distribute verbatim copies
 of this license document, but changing it is not allowed.

@ -278,3 +278,61 @@ PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
 POSSIBILITY OF SUCH DAMAGES.

                     END OF TERMS AND CONDITIONS
+
+            How to Apply These Terms to Your New Programs
+
+  If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+  To do so, attach the following notices to the program.  It is safest
+to attach them to the start of each source file to most effectively
+convey the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+    <one line to give the program's name and a brief idea of what it does.>
+    Copyright (C) <year>  <name of author>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, see <https://www.gnu.org/licenses/>.
+
+Also add information on how to contact you by electronic and paper mail.
+
+If the program is interactive, make it output a short notice like this
+when it starts in an interactive mode:
+
+    Gnomovision version 69, Copyright (C) year name of author
+    Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+    This is free software, and you are welcome to redistribute it
+    under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License.  Of course, the commands you use may
+be called something other than `show w' and `show c'; they could even be
+mouse-clicks or menu items--whatever suits your program.
+
+You should also get your employer (if you work as a programmer) or your
+school, if any, to sign a "copyright disclaimer" for the program, if
+necessary.  Here is a sample; alter the names:
+
+  Yoyodyne, Inc., hereby disclaims all copyright interest in the program
+  `Gnomovision' (which makes passes at compilers) written by James Hacker.
+
+  <signature of Moe Ghoul>, 1 April 1989
+  Moe Ghoul, President of Vice
+
+This General Public License does not permit incorporating your program into
+proprietary programs.  If your program is a subroutine library, you may
+consider it more useful to permit linking proprietary applications with the
+library.  If this is what you want to do, use the GNU Lesser General
+Public License instead of this License.
--- a/README.md
+++ b/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 -
-  <https://github.com/usgs-coupled/iphreeqc> -
-  <https://git.gfz-potsdam.de/naaice/iphreeqc>
- **tug** - <https://git.gfz-potsdam.de/naaice/tug>
+- **litephreeqc**: IPhreeqc
+   (<https://github.com/usgs-coupled/iphreeqc>) with patches from
+   GFZ/UP: <https://git.gfz.de/naaice/litephreeqc>
+- **tug** - <https://git.gfz.de/naaice/tug>

 ## Installation

@ -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
--- a/R_lib/init_r_lib.R
+++ b/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
--- a/R_lib/kin_r_library.R
+++ b/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"))
+## }
--- a/bench/CMakeLists.txt
+++ b/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}
--- a/bench/barite/barite_200.R
+++ b/bench/barite/barite_200.R
@ -35,15 +35,15 @@ diffusion_setup <- list(
 )

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

 chemistry_setup <- list(
--- a/bench/barite/barite_50ai.R
+++ b/bench/barite/barite_50ai.R
@ -11,9 +11,9 @@ grid_def <- matrix(2, nrow = rows, ncol = cols)
 grid_setup <- list(
  pqc_in_file = "./barite.pqi",
  pqc_db_file = "./db_barite.dat", ## Path to the database file for Phreeqc
-  grid_def = grid_def,             ## Definition of the grid, containing IDs according to the Phreeqc input script
-  grid_size = c(s_rows, s_cols),   ## Size of the grid in meters
-  constant_cells = c()             ## IDs of cells with constant concentration
+  grid_def = grid_def, ## Definition of the grid, containing IDs according to the Phreeqc input script
+  grid_size = c(s_rows, s_cols), ## Size of the grid in meters
+  constant_cells = c() ## IDs of cells with constant concentration
 )

 bound_length <- 2
@ -36,15 +36,15 @@ diffusion_setup <- list(
 )

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

 chemistry_setup <- list(
--- a/bench/dolo/dolo_interp.R
+++ b/bench/dolo/dolo_interp.R
@ -41,36 +41,36 @@ diffusion_setup <- list(
 )

 check_sign_cal_dol_dht <- function(old, new) {
-    if ((old["Calcite"] == 0) != (new["Calcite"] == 0)) {
-        return(TRUE)
-    }
-    if ((old["Dolomite"] == 0) != (new["Dolomite"] == 0)) {
-        return(TRUE)
-    }
+    # if ((old["Calcite"] == 0) != (new["Calcite"] == 0)) {
+    #     return(TRUE)
+    # }
+    # if ((old["Dolomite"] == 0) != (new["Dolomite"] == 0)) {
+    #     return(TRUE)
+    # }
    return(FALSE)
 }

-fuzz_input_dht_keys <- function(input) {
-    dht_species <- c(
-        "H" = 3,
-        "O" = 3,
-        "Charge" = 3,
-        "C(4)" = 6,
-        "Ca" = 6,
-        "Cl" = 3,
-        "Mg" = 5,
-        "Calcite" = 4,
-        "Dolomite" = 4
-    )
-    return(input[names(dht_species)])
-}
+# fuzz_input_dht_keys <- function(input) {
+#     dht_species <- c(
+#         "H" = 3,
+#         "O" = 3,
+#         "Charge" = 3,
+#         "C" = 6,
+#         "Ca" = 6,
+#         "Cl" = 3,
+#         "Mg" = 5,
+#         "Calcite" = 4,
+#         "Dolomite" = 4
+#     )
+#     return(input[names(dht_species)])
+# }

 check_sign_cal_dol_interp <- function(to_interp, data_set) {
    dht_species <- c(
        "H" = 3,
        "O" = 3,
        "Charge" = 3,
-        "C(4)" = 6,
+        "C" = 6,
        "Ca" = 6,
        "Cl" = 3,
        "Mg" = 5,
@ -95,7 +95,7 @@ check_neg_cal_dol <- function(result) {
 # significant digits)

 pht_species <- c(
-    "C(4)" = 3,
+    "C" = 3,
    "Ca" = 3,
    "Mg" = 2,
    "Calcite" = 2,
@ -107,7 +107,7 @@ chemistry_setup <- list(
        "H" = 3,
        "O" = 3,
        "Charge" = 3,
-        "C(4)" = 6,
+        "C" = 6,
        "Ca" = 6,
        "Cl" = 3,
        "Mg" = 5,
@ -117,7 +117,7 @@ chemistry_setup <- list(
    pht_species = pht_species,
    hooks = list(
        dht_fill = check_sign_cal_dol_dht,
-        dht_fuzz = fuzz_input_dht_keys,
+        # dht_fuzz = fuzz_input_dht_keys,
        interp_pre = check_sign_cal_dol_interp,
        interp_post = check_neg_cal_dol
    )
--- a/bench/fgcs/20241211_README.pdf
+++ b/bench/fgcs/20241211_README.pdf
--- a/bench/fgcs/20241211_README.tex
+++ b/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:
--- a/bench/fgcs/EvalFGCS.R
+++ b/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
--- a/bench/fgcs/README.org
+++ b/bench/fgcs/README.org
@ -0,0 +1,2 @@
+* Refer to the LaTeX file (and pdf) for more information
+
--- a/bench/fgcs/barite_fgcs_2.R
+++ b/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
+)
--- a/bench/fgcs/barite_fgcs_2.pqi
+++ b/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
--- a/bench/fgcs/barite_fgcs_2_rt.R
+++ b/bench/fgcs/barite_fgcs_2_rt.R
@ -0,0 +1,7 @@
+iterations <- 200
+dt <- 1000
+
+list(
+    timesteps = rep(dt, iterations),
+    store_result = TRUE
+)
--- a/docs/POET.drawio
+++ b/docs/POET.drawio
--- a/docs/POET_scheme.svg
+++ b/docs/POET_scheme.svg
--- a/ext/iphreeqc
+++ b/ext/iphreeqc
@ -1 +0,0 @@
-Subproject commit e6e5e0d5156c093241a53e6ce074ef346d64ae26
--- a/ext/litephreeqc
+++ b/ext/litephreeqc
@ -0,0 +1 @@
+Subproject commit 953c752431d2b2758268083f407f943843efc7ad
--- a/ext/tug
+++ b/ext/tug
@ -1 +1 @@
-Subproject commit 449647010ab9cdf9e405139f360424a2b21ab3ab
+Subproject commit 9c4aeee410c71d064f7567143d4f8d6451ade75a
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@ -1,7 +1,4 @@
-add_library(POETLib)
-
-target_sources(POETLib 
-  PRIVATE
+add_library(POETLib
  Init/InitialList.cpp 
  Init/GridInit.cpp 
  Init/DiffusionInit.cpp
@ -18,6 +15,15 @@ target_sources(POETLib
  Chemistry/SurrogateModels/ProximityHashTable.cpp
 )

+set(POET_TUG_APPROACH "Implicit" CACHE STRING "tug numerical approach to use")
+set_property(CACHE POET_TUG_APPROACH PROPERTY STRINGS "Implicit" "Explicit")
+
+if (POET_TUG_APPROACH STREQUAL "Implicit")
+  target_compile_definitions(POETLib PRIVATE POET_TUG_BTCS)
+elseif (POET_TUG_APPROACH STREQUAL "Explicit")
+  target_compile_definitions(POETLib PRIVATE POET_TUG_FTCS)
+endif()
+
 target_include_directories(POETLib PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}")
 target_link_libraries(
  POETLib 
@ -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)
--- a/src/Chemistry/ChemistryDefs.hpp
+++ b/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) {
-    input.resize(size);
-    output.resize(size);
-    mapping.resize(size, CHEM_PQC);
-  }
+  WorkPackage(std::size_t _size)
+      : size(_size), input(size), output(size), mapping(size, CHEM_PQC) {}
 };
 } // namespace poet
--- a/src/Chemistry/ChemistryModule.cpp
+++ b/src/Chemistry/ChemistryModule.cpp
@ -1,11 +1,14 @@
 #include "ChemistryModule.hpp"

 #include "PhreeqcEngine.hpp"
+#include "PhreeqcMatrix.hpp"
+#include "PhreeqcRunner.hpp"
 #include "SurrogateModels/DHT_Wrapper.hpp"
 #include "SurrogateModels/Interpolation.hpp"

 #include <algorithm>
 #include <cassert>
+#include <cmath>
 #include <cstdint>
 #include <memory>
 #include <mpi.h>
@ -63,7 +66,8 @@ inverseDistanceWeighting(const std::vector<std::int32_t> &to_calc,
      distance += std::pow(
          rescaled[key_comp_i][point_i] - rescaled[key_comp_i][data_set_n], 2);
    }
-    weights[point_i] = 1 / std::sqrt(distance);
+
+    weights[point_i] = distance != 0 ? 1 / std::sqrt(distance) : 0;
    assert(!std::isnan(weights[point_i]));
    inv_sum += weights[point_i];
  }
@ -94,63 +98,9 @@ inverseDistanceWeighting(const std::vector<std::int32_t> &to_calc,
    key_delta /= inv_sum;

    results[output_comp_i] = from[output_comp_i] + key_delta;
+    assert(!std::isnan(results[output_comp_i]));
  }

-  // if (!has_h) {
-  //   double new_val = 0;
-  //   for (int j = 0; j < data_set_n; j++) {
-  //     new_val += weights[j] * output[j][0];
-  //   }
-  //   results[0] = new_val / inv_sum;
-  // }
-
-  // if (!has_h) {
-  //   double new_val = 0;
-  //   for (int j = 0; j < data_set_n; j++) {
-  //     new_val += weights[j] * output[j][1];
-  //   }
-  //   results[1] = new_val / inv_sum;
-  // }
-
-  // for (std::uint32_t i = 0; i < to_calc.size(); i++) {
-  //   const std::uint32_t interp_i = to_calc[i];
-
-  //   // rescale input between 0 and 1
-  //   for (int j = 0; j < input.size(); j++) {
-  //     buffer[j] = input[j].at(i);
-  //   }
-
-  //   buffer[buffer_size - 1] = from[interp_i];
-
-  //   const double min = *std::min_element(buffer, buffer + buffer_size);
-  //   const double max = *std::max_element(buffer, buffer + buffer_size);
-
-  //   for (int j = 0; j < input.size(); j++) {
-  //     buffer[j] = ((max - min) != 0 ? (buffer[j] - min) / (max - min) : 1);
-  //   }
-  //   from_rescaled =
-  //       ((max - min) != 0 ? (from[interp_i] - min) / (max - min) : 0);
-
-  //   double inv_sum = 0;
-
-  //   // calculate distances for each point
-  //   for (int i = 0; i < input.size(); i++) {
-  //     const double distance = std::pow(buffer[i] - from_rescaled, 2);
-
-  //     buffer[i] = distance > 0 ? (1 / std::sqrt(distance)) : 0;
-  //     inv_sum += buffer[i];
-  //   }
-  //   // calculate new values
-  //   double new_val = 0;
-  //   for (int i = 0; i < output.size(); i++) {
-  //     new_val += buffer[i] * output[i][interp_i];
-  //   }
-  //   results[interp_i] = new_val / inv_sum;
-  //   if (std::isnan(results[interp_i])) {
-  //     std::cout << "nan with new_val = " << output[0][i] << std::endl;
-  //   }
-  // }
-
  return results;
 }

@ -167,16 +117,12 @@ poet::ChemistryModule::ChemistryModule(
  this->n_cells = chem_params.total_grid_cells;

  if (!is_master) {
-    for (const auto &[pqc_id, pqc_config] : chem_params.pqc_config) {
-      this->phreeqc_instances[pqc_id] =
-          std::make_unique<PhreeqcEngine>(pqc_config);
-    }
-    // for (std::size_t i = 0; i < chem_params.pqc_ids.size(); i++) {
-    //   this->phreeqc_instances[chem_params.pqc_ids[i]] =
-    //       std::make_unique<PhreeqcWrapper>(
-    //           chem_params.database, chem_params.pqc_scripts[i],
-    //           chem_params.pqc_sol_order, chem_params.field_header, wp_size_);
-    // }
+    PhreeqcMatrix pqc_mat =
+        PhreeqcMatrix(chem_params.database, chem_params.pqc_script,
+                      chem_params.with_h0_o0, chem_params.with_redox);
+
+    this->pqc_runner =
+        std::make_unique<PhreeqcRunner>(pqc_mat.subset(chem_params.pqc_ids));
  }
 }

@ -187,11 +133,10 @@ poet::ChemistryModule::~ChemistryModule() {
 }

 void poet::ChemistryModule::initializeDHT(
-    uint32_t size_mb, const NamedVector<std::uint32_t> &key_species) {
+    uint32_t size_mb, const NamedVector<std::uint32_t> &key_species,
+    bool has_het_ids) {
  constexpr uint32_t MB_FACTOR = 1E6;

-  this->dht_enabled = true;
-
  MPI_Comm dht_comm;

  if (this->is_master) {
@ -221,7 +166,7 @@ void poet::ChemistryModule::initializeDHT(

    this->dht = new DHT_Wrapper(dht_comm, dht_size, map_copy, key_indices,
                                this->prop_names, params.hooks,
-                                this->prop_count, interp_enabled);
+                                this->prop_count, interp_enabled, has_het_ids);
    this->dht->setBaseTotals(base_totals.at(0), base_totals.at(1));
  }
 }
@ -312,9 +257,10 @@ void poet::ChemistryModule::initializeInterp(
      map_copy = this->dht->getKeySpecies();
      for (auto i = 0; i < map_copy.size(); i++) {
        const std::uint32_t signif =
-            static_cast<std::uint32_t>(map_copy[i]) - (map_copy[i] > InterpolationModule::COARSE_DIFF
-                               ? InterpolationModule::COARSE_DIFF
-                               : 0);
+            static_cast<std::uint32_t>(map_copy[i]) -
+            (map_copy[i] > InterpolationModule::COARSE_DIFF
+                 ? InterpolationModule::COARSE_DIFF
+                 : 0);
        map_copy[i] = signif;
      }
    }
@ -371,7 +317,8 @@ void poet::ChemistryModule::unshuffleField(const std::vector<double> &in_buffer,
    }
  }
 }
-  
-void poet::ChemistryModule::set_ai_surrogate_validity_vector(std::vector<int> r_vector) {
+
+void poet::ChemistryModule::set_ai_surrogate_validity_vector(
+    std::vector<int> r_vector) {
  this->ai_surrogate_validity_vector = r_vector;
 }
--- a/src/Chemistry/ChemistryModule.hpp
+++ b/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) {
@ -222,8 +234,8 @@ public:
  };

  /**
-  *  **Master only** Set the ai surrogate validity vector from R
-  */
+   *  **Master only** Set the ai surrogate validity vector from R
+   */
  void set_ai_surrogate_validity_vector(std::vector<int> r_vector);

  std::vector<uint32_t> GetWorkerInterpolationCalls() const;
@ -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

--- a/src/Chemistry/SurrogateModels/DHT.c
+++ b/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;
--- a/src/Chemistry/SurrogateModels/DHT.h
+++ b/src/Chemistry/SurrogateModels/DHT.h
@ -117,6 +117,9 @@ typedef struct {
  unsigned int index_count;

  int (*accumulate_callback)(int, void *, int, void *);
+
+  size_t sum_idx;
+  size_t cnt_idx;
 #ifdef DHT_STATISTICS
  /** Detailed statistics of the usage of the DHT. */
  DHT_stats *stats;
@ -125,10 +128,11 @@ typedef struct {

 extern void DHT_set_accumulate_callback(DHT *table,
                                        int (*callback_func)(int, void *, int,
-                                                              void *));
+                                                             void *));

 extern int DHT_write_accumulate(DHT *table, const void *key, int send_size,
-                                void *data, uint32_t *proc, uint32_t *index, int *callback_ret);
+                                void *data, uint32_t *proc, uint32_t *index,
+                                int *callback_ret);

 /**
 * @brief Create a DHT.
@ -284,44 +288,8 @@ extern int DHT_free(DHT *table, int *eviction_counter, int *readerror_counter);
 */
 extern int DHT_print_statistics(DHT *table);

-/**
- * @brief Determine destination rank and index.
- *
- * This is done by looping over all possbile indices. First of all, set a
- * temporary index to zero and copy count of bytes for each index into the
- * memory area of the temporary index. After that the current index is
- * calculated by the temporary index modulo the table size. The destination rank
- * of the process is simply determined by hash modulo the communicator size.
- *
- * @param hash Calculated 64 bit hash.
- * @param comm_size Communicator size.
- * @param table_size Count of buckets per process.
- * @param dest_rank Reference to the destination rank variable.
- * @param index Pointer to the array index.
- * @param index_count Count of possible indeces.
- */
-static void determine_dest(uint64_t hash, int comm_size,
-                           unsigned int table_size, unsigned int *dest_rank,
-                           uint64_t *index, unsigned int index_count);
+extern float DHT_get_used_idx_factor(DHT *table, int with_reset);

-/**
- * @brief Set the occupied flag.
- *
- * This will set the first bit of a bucket to 1.
- *
- * @param flag_byte First byte of a bucket.
- */
-static void set_flag(char *flag_byte);
-
-/**
- * @brief Get the occupied flag.
- *
- * This function determines whether the occupied flag of a bucket was set or
- * not.
- *
- * @param flag_byte First byte of a bucket.
- * @return int Returns 1 for true or 0 for false.
- */
-static int read_flag(char flag_byte);
+extern int DHT_flush(DHT *table);

 #endif /* DHT_H */
--- a/src/Chemistry/SurrogateModels/DHT_Wrapper.cpp
+++ b/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,42 +129,47 @@ void DHT_Wrapper::fillDHT(const WorkPackage &work_package) {
  dht_results.filledDHT = std::vector<bool>(length, false);
  for (int i = 0; i < length; i++) {
    // If true grid cell was simulated, needs to be inserted into dht
-    if (work_package.mapping[i] == CHEM_PQC) {
-
-      // check if calcite or dolomite is absent and present, resp.n and vice
-      // versa in input/output. If this is the case -> Do not write to DHT!
-      // HACK: hardcoded, should be fixed!
-      if (hooks.dht_fill.isValid()) {
-        NamedVector<double> old_values(output_names, work_package.input[i]);
-        NamedVector<double> new_values(output_names, work_package.output[i]);
-
-        if (hooks.dht_fill(old_values, new_values)) {
-          continue;
-        }
-      }
-
-      uint32_t proc, index;
-      auto &key = dht_results.keys[i];
-      const auto data =
-          (with_interp ? outputToInputAndRates(work_package.input[i],
-                                               work_package.output[i])
-                       : work_package.output[i]);
-      // void *data = (void *)&(work_package[i * this->data_count]);
-      // fuzz data (round, logarithm etc.)
-
-      // insert simulated data with fuzzed key into DHT
-      int res = DHT_write(this->dht_object, key.data(),
-                          const_cast<double *>(data.data()), &proc, &index);
-
-      dht_results.locations[i] = {proc, index};
-
-      // if data was successfully written ...
-      if ((res != DHT_SUCCESS) && (res == DHT_WRITE_SUCCESS_WITH_EVICTION)) {
-        dht_evictions++;
-      }
-
-      dht_results.filledDHT[i] = true;
+    if (work_package.mapping[i] != CHEM_PQC) {
+      continue;
    }
+
+    if (work_package.input[i][0] != 2) {
+      continue;
+    }
+
+    // check if calcite or dolomite is absent and present, resp.n and vice
+    // versa in input/output. If this is the case -> Do not write to DHT!
+    // HACK: hardcoded, should be fixed!
+    if (hooks.dht_fill.isValid()) {
+      NamedVector<double> old_values(output_names, work_package.input[i]);
+      NamedVector<double> new_values(output_names, work_package.output[i]);
+
+      if (Rcpp::as<bool>(hooks.dht_fill(old_values, new_values))) {
+        continue;
+      }
+    }
+
+    uint32_t proc, index;
+    auto &key = dht_results.keys[i];
+    const auto data =
+        (with_interp ? outputToInputAndRates(work_package.input[i],
+                                             work_package.output[i])
+                     : work_package.output[i]);
+    // void *data = (void *)&(work_package[i * this->data_count]);
+    // fuzz data (round, logarithm etc.)
+
+    // insert simulated data with fuzzed key into DHT
+    int res = DHT_write(this->dht_object, key.data(),
+                        const_cast<double *>(data.data()), &proc, &index);
+
+    dht_results.locations[i] = {proc, index};
+
+    // if data was successfully written ...
+    if ((res != DHT_SUCCESS) && (res == DHT_WRITE_SUCCESS_WITH_EVICTION)) {
+      dht_evictions++;
+    }
+
+    dht_results.filledDHT[i] = true;
  }
 }

@ -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;
 }
--- a/src/Chemistry/SurrogateModels/DHT_Wrapper.hpp
+++ b/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

--- a/src/Chemistry/SurrogateModels/Interpolation.hpp
+++ b/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

--- a/src/Chemistry/SurrogateModels/InterpolationModule.cpp
+++ b/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;
      }
--- a/src/Chemistry/SurrogateModels/LookupKey.hpp
+++ b/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> {
--- a/src/Chemistry/SurrogateModels/Rounding.hpp
+++ b/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)))) -
--- a/src/Chemistry/WorkerFunctions.cpp
+++ b/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);
    }
+  }
+  this->pqc_runner->run(inout_chem, dTimestep, to_ignore);

-    auto curr_input = work_package.input[wp_id];
-    const auto pqc_id = static_cast<int>(curr_input[0]);
-
-    auto &phreeqc_instance = this->phreeqc_instances[pqc_id];
-    work_package.output[wp_id] = work_package.input[wp_id];
-
-    curr_input.erase(std::remove_if(curr_input.begin(), curr_input.end(),
-                                    [](double d) { return std::isnan(d); }),
-                     curr_input.end());
-
-    phreeqc_instance->runCell(curr_input, dTimestep);
-
-    std::size_t output_index = 0;
-    for (std::size_t i = 0; i < work_package.output[wp_id].size(); i++) {
-      if (!std::isnan(work_package.output[wp_id][i])) {
-        work_package.output[wp_id][i] = curr_input[output_index++];
-      }
+  for (std::size_t wp_id = 0; wp_id < work_package.size; wp_id++) {
+    if (work_package.mapping[wp_id] == CHEM_PQC) {
+      work_package.output[wp_id] = inout_chem[wp_id];
    }
  }
 }
--- a/src/Init/ChemistryInit.cpp
+++ b/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++) {
-    POETInitCell cell = {
-        pqc_solutions,
-        pqc_solution_primaries,
-        Rcpp::as<std::vector<std::string>>(pqc_exchanger[i]),
-        Rcpp::as<std::vector<std::string>>(pqc_kinetics[i]),
-        Rcpp::as<std::vector<std::string>>(pqc_equilibrium[i]),
-        Rcpp::as<std::vector<std::string>>(pqc_surface_comps[i]),
-        Rcpp::as<std::vector<std::string>>(pqc_surface_charges[i])};
-
-    chem_init.pqc_config[pqc_ids[i]] = {database, pqc_scripts[i], cell};
-  }
+  // for (std::size_t i = 0; i < pqc_ids.size(); ++i) {
+  //   chem_init.pqc_input[pqc_ids[i]] = pqc_scripts[i];
+  // }

  // chem_init.pqc_sol_order = pqc_solutions;

--- a/src/Init/DiffusionInit.cpp
+++ b/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>
-extend_transport_names(std::unique_ptr<PhreeqcInit> &phreeqc,
-                       const Rcpp::List &boundaries_list,
-                       const Rcpp::List &inner_boundaries,
-                       const std::vector<std::string> &old_trans_names) {
+// static std::vector<std::string>
+// extend_transport_names(std::unique_ptr<PhreeqcInit> &phreeqc,
+//                        const Rcpp::List &boundaries_list,
+//                        const Rcpp::List &inner_boundaries,
+//                        const std::vector<std::string> &old_trans_names) {

-  std::vector<std::string> transport_names = old_trans_names;
-  std::set<int> constant_pqc_ids;
+//   std::vector<std::string> transport_names = old_trans_names;
+//   std::set<int> constant_pqc_ids;

-  for (const auto &side : tug_side_mapping) {
-    if (!boundaries_list.containsElementNamed(side.second.c_str())) {
-      continue;
-    }
+//   for (const auto &side : tug_side_mapping) {
+//     if (!boundaries_list.containsElementNamed(side.second.c_str())) {
+//       continue;
+//     }

-    Rcpp::List mapping = boundaries_list[side.second];
+//     Rcpp::List mapping = boundaries_list[side.second];

-    const Rcpp::NumericVector cells = mapping["cell"];
-    const Rcpp::NumericVector values = mapping["sol_id"];
-    const Rcpp::CharacterVector type_str = mapping["type"];
+//     const Rcpp::NumericVector cells = mapping["cell"];
+//     const Rcpp::NumericVector values = mapping["sol_id"];
+//     const Rcpp::CharacterVector type_str = mapping["type"];

-    if (cells.size() != values.size()) {
-      throw std::runtime_error("Boundary [" + side.second +
-                               "] cells and values are not the same "
-                               "length");
-    }
+//     if (cells.size() != values.size()) {
+//       throw std::runtime_error("Boundary [" + side.second +
+//                                "] cells and values are not the same "
+//                                "length");
+//     }

-    for (auto i = 0; i < cells.size(); i++) {
-      if (type_str[i] == "constant") {
-        constant_pqc_ids.insert(static_cast<std::size_t>(values[i]));
-      }
-    }
-  }
+//     for (auto i = 0; i < cells.size(); i++) {
+//       if (type_str[i] == "constant") {
+//         constant_pqc_ids.insert(static_cast<std::size_t>(values[i]));
+//       }
+//     }
+//   }

-  if (inner_boundaries.size() > 0) {
-    const Rcpp::NumericVector values = inner_boundaries["sol_id"];
-    for (auto i = 0; i < values.size(); i++) {
-      constant_pqc_ids.insert(static_cast<std::size_t>(values[i]));
-    }
-  }
+//   if (inner_boundaries.size() > 0) {
+//     const Rcpp::NumericVector values = inner_boundaries["sol_id"];
+//     for (auto i = 0; i < values.size(); i++) {
+//       constant_pqc_ids.insert(static_cast<std::size_t>(values[i]));
+//     }
+//   }

-  if (!constant_pqc_ids.empty()) {
-    constexpr std::size_t keep_h_o_charge = 3;
+//   if (!constant_pqc_ids.empty()) {
+//     constexpr std::size_t keep_h_o_charge = 3;

-    for (const auto &pqc_id : constant_pqc_ids) {
-      const auto solution_names = phreeqc->getSolutionNames(pqc_id);
+//     for (const auto &pqc_id : constant_pqc_ids) {
+//       const auto solution_names = phreeqc->getSolutionNames(pqc_id);

-      // add those strings which are not already in the transport_names
-      for (const auto &name : solution_names) {
-        if (std::find(transport_names.begin(), transport_names.end(), name) ==
-            transport_names.end()) {
-          auto position =
-              std::lower_bound(transport_names.begin() + keep_h_o_charge,
-                               transport_names.end(), name);
+//       // add those strings which are not already in the transport_names
+//       for (const auto &name : solution_names) {
+//         if (std::find(transport_names.begin(), transport_names.end(), name)
+//         ==
+//             transport_names.end()) {
+//           auto position =
+//               std::lower_bound(transport_names.begin() + keep_h_o_charge,
+//                                transport_names.end(), name);

-          transport_names.insert(position, name);
-        }
-      }
-    }
-  }
+//           transport_names.insert(position, name);
+//         }
+//       }
+//     }
+//   }

-  return transport_names;
-}
+//   return transport_names;
+// }

-static Rcpp::List
-extend_initial_grid(const Rcpp::List &initial_grid,
-                    const std::vector<std::string> &transport_names) {
-  // std::vector<std::string> names_to_add(transport_names.begin() + old_size,
-  //                                       transport_names.end());
+// static Rcpp::List
+// extend_initial_grid(const Rcpp::List &initial_grid,
+//                     const std::vector<std::string> &transport_names) {
+//   // std::vector<std::string> names_to_add(transport_names.begin() +
+//   old_size,
+//   //                                       transport_names.end());

-  Rcpp::Function extend_grid_R("add_missing_transport_species");
+//   Rcpp::Function extend_grid_R("add_missing_transport_species");

-  return extend_grid_R(initial_grid, Rcpp::wrap(transport_names));
-}
+//   return extend_grid_R(initial_grid, Rcpp::wrap(transport_names));
+// }

 std::pair<Rcpp::List, Rcpp::List>
 InitialList::resolveBoundaries(const Rcpp::List &boundaries_list,
-                               const Rcpp::List &inner_boundaries) {
+                               const Rcpp::List &inner_boundaries,
+                               const PhreeqcMatrix &phreeqc_mat) {
  Rcpp::List bound_list;
  Rcpp::List inner_bound;
  Rcpp::Function resolve_R("resolve_pqc_bound");

-  const std::size_t old_transport_size = this->transport_names.size();
+  // const std::size_t old_transport_size = this->transport_names.size();

-  this->transport_names = extend_transport_names(
-      this->phreeqc, boundaries_list, inner_boundaries, this->transport_names);
+  // this->transport_names = extend_transport_names(
+  //     this->phreeqc, boundaries_list, inner_boundaries,
+  //     this->transport_names);

-  const std::size_t new_transport_size = this->transport_names.size();
+  // const std::size_t new_transport_size = this->transport_names.size();

-  if (old_transport_size != new_transport_size) {
-    this->initial_grid =
-        extend_initial_grid(this->initial_grid, this->transport_names);
-  }
+  // if (old_transport_size != new_transport_size) {
+  //   this->initial_grid =
+  //       extend_initial_grid(this->initial_grid, this->transport_names);
+  // }

  for (const auto &species_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>(
-                resolve_R(this->phreeqc_mat, Rcpp::wrap(species_name), values[i]));
+            c_value[c_id] =
+                static_cast<TugType>(phreeqc_mat(values[i], species_name));
+            // c_value[c_id] = Rcpp::as<TugType>(
+            //     resolve_R(this->phreeqc_mat, Rcpp::wrap(species_name),
+            //     values[i]));
          } else {
            throw std::runtime_error("Unknown boundary type");
          }
@ -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(
-            this->phreeqc_mat, Rcpp::wrap(species_name), inner_pqc_id_vec[i])));
+        c_value.push_back(static_cast<TugType>(
+            phreeqc_mat(inner_pqc_id_vec[i], species_name)));
+        // c_value.push_back(Rcpp::as<TugType>(resolve_R(
+        //     this->phreeqc_mat, Rcpp::wrap(species_name),
+        //     inner_pqc_id_vec[i])));
      }

      inner_bound[species_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;

--- a/src/Init/GridInit.cpp
+++ b/src/Init/GridInit.cpp
@ -1,13 +1,19 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+/*
+ * GridInit.cpp - Implementation of grid initialization
+ * Copyright (C) 2024-2025 Max Luebke (University of Potsdam)
+ */
+
 #include "InitialList.hpp"

-#include "PhreeqcInit.hpp"
+#include "PhreeqcMatrix.hpp"

 #include <RInside.h>
 #include <Rcpp/Function.h>
+#include <Rcpp/vector/Matrix.h>
 #include <Rcpp/vector/instantiation.h>
-#include <cstdint>
+#include <fstream>
 #include <map>
-#include <memory>
 #include <regex>
 #include <sstream>
 #include <string>
@ -15,38 +21,6 @@

 namespace poet {

-static Rcpp::NumericMatrix
-pqcScriptToGrid(std::unique_ptr<PhreeqcInit> &phreeqc, RInside &R) {
-  PhreeqcInit::PhreeqcMat phreeqc_mat = phreeqc->getPhreeqcMat();
-
-  // add "id" to the front of the column names
-
-  const std::size_t ncols = phreeqc_mat.names.size();
-  const std::size_t nrows = phreeqc_mat.values.size();
-
-  phreeqc_mat.names.insert(phreeqc_mat.names.begin(), "ID");
-
-  Rcpp::NumericMatrix mat(nrows, ncols + 1);
-
-  for (std::size_t i = 0; i < nrows; i++) {
-    mat(i, 0) = phreeqc_mat.ids[i];
-    for (std::size_t j = 0; j < ncols; ++j) {
-      mat(i, j + 1) = phreeqc_mat.values[i][j];
-    }
-  }
-
-  Rcpp::colnames(mat) = Rcpp::wrap(phreeqc_mat.names);
-
-  return mat;
-}
-
-static inline Rcpp::List matToGrid(RInside &R, const Rcpp::NumericMatrix &mat,
-                                   const Rcpp::NumericMatrix &grid) {
-  Rcpp::Function pqc_to_grid_R("pqc_to_grid");
-
-  return pqc_to_grid_R(mat, grid);
-}
-
 static inline std::map<int, std::string>
 replaceRawKeywordIDs(std::map<int, std::string> raws) {
  for (auto &raw : raws) {
@ -60,25 +34,6 @@ replaceRawKeywordIDs(std::map<int, std::string> raws) {
  return raws;
 }

-static inline uint32_t getSolutionCount(std::unique_ptr<PhreeqcInit> &phreeqc,
-                                        const Rcpp::List &initial_grid) {
-  PhreeqcInit::ModulesArray mod_array;
-  Rcpp::Function unique_R("unique");
-
-  std::vector<int> row_ids =
-      Rcpp::as<std::vector<int>>(unique_R(initial_grid["ID"]));
-
-  // std::vector<std::uint32_t> sizes_vec(sizes.begin(), sizes.end());
-
-  // Rcpp::Function modify_sizes("modify_module_sizes");
-  // sizes_vec = Rcpp::as<std::vector<std::uint32_t>>(
-  //     modify_sizes(sizes_vec, phreeqc_mat, initial_grid));
-
-  // std::copy(sizes_vec.begin(), sizes_vec.end(), sizes.begin());
-
-  return phreeqc->getModuleSizes(row_ids)[POET_SOL];
-}
-
 static std::string readFile(const std::string &path) {
  std::string string_rpath(PATH_MAX, '\0');

@ -98,10 +53,26 @@ static std::string readFile(const std::string &path) {
  return buffer.str();
 }

-void InitialList::prepareGrid(const Rcpp::List &grid_input) {
-  // parse input values
-  Rcpp::Function unique_R("unique");
+static Rcpp::List expandGrid(const PhreeqcMatrix &pqc_mat,
+                             const std::vector<int> unique_ids,
+                             const Rcpp::IntegerMatrix &grid_def) {

+  PhreeqcMatrix subset_pqc_mat = pqc_mat.subset(unique_ids);
+
+  PhreeqcMatrix::STLExport phreeqc_mat =
+      subset_pqc_mat.get(PhreeqcMatrix::VectorExportType::COLUMN_MAJOR);
+
+  const std::size_t ncols = phreeqc_mat.names.size();
+  const std::size_t nrows = phreeqc_mat.values.size() / ncols;
+
+  Rcpp::NumericMatrix pqc_mat_R(nrows, ncols, phreeqc_mat.values.begin());
+  Rcpp::colnames(pqc_mat_R) = Rcpp::wrap(phreeqc_mat.names);
+
+  return Rcpp::Function("pqc_to_grid")(pqc_mat_R, grid_def);
+}
+
+PhreeqcMatrix InitialList::prepareGrid(const Rcpp::List &grid_input) {
+  // parse input values
  std::string script;
  std::string database;

@ -126,8 +97,9 @@ void InitialList::prepareGrid(const Rcpp::List &grid_input) {
  }

  this->database = database;
+  this->pqc_script = script;

-  Rcpp::NumericMatrix grid_def =
+  Rcpp::IntegerMatrix grid_def =
      grid_input[GRID_MEMBER_STR(GridMembers::GRID_DEF)];
  Rcpp::NumericVector grid_size =
      grid_input[GRID_MEMBER_STR(GridMembers::GRID_SIZE)];
@ -156,35 +128,77 @@ void InitialList::prepareGrid(const Rcpp::List &grid_input) {
    throw std::runtime_error("Grid size must be positive.");
  }

-  this->phreeqc = std::make_unique<PhreeqcInit>(database, script);
+  bool with_redox =
+      grid_input.containsElementNamed(
+          GRID_MEMBER_STR(GridMembers::PQC_WITH_REDOX))
+          ? Rcpp::as<bool>(
+                grid_input[GRID_MEMBER_STR(GridMembers::PQC_WITH_REDOX)])
+          : false;

-  this->phreeqc_mat = pqcScriptToGrid(phreeqc, R);
-  this->initial_grid = matToGrid(R, this->phreeqc_mat, grid_def);
+  bool with_h0_o0 =
+      grid_input.containsElementNamed(
+          GRID_MEMBER_STR(GridMembers::PQC_WITH_H0_O0))
+          ? Rcpp::as<bool>(
+                grid_input[GRID_MEMBER_STR(GridMembers::PQC_WITH_H0_O0)])
+          : false;

-  const uint32_t solution_count = getSolutionCount(phreeqc, this->initial_grid);
-
-  std::vector<std::string> colnames =
-      Rcpp::as<std::vector<std::string>>(this->initial_grid.names());
-
-  this->transport_names = std::vector<std::string>(
-      colnames.begin() + 1,
-      colnames.begin() + 1 + solution_count); // skip ID
-
-  std::map<int, std::string> pqc_raw_dumps;
-
-  pqc_raw_dumps = replaceRawKeywordIDs(phreeqc->raw_dumps());
-
-  this->pqc_ids =
-      Rcpp::as<std::vector<int>>(unique_R(this->initial_grid["ID"]));
-
-  for (const auto &id : this->pqc_ids) {
-    this->pqc_scripts.push_back(pqc_raw_dumps[id]);
-    this->pqc_exchanger.push_back(phreeqc->getExchanger(id));
-    this->pqc_kinetics.push_back(phreeqc->getKineticsNames(id));
-    this->pqc_equilibrium.push_back(phreeqc->getEquilibriumNames(id));
-    this->pqc_surface_comps.push_back(phreeqc->getSurfaceCompNames(id));
-    this->pqc_surface_charges.push_back(phreeqc->getSurfaceChargeNames(id));
+  if (with_h0_o0 && !with_redox) {
+    throw std::runtime_error(
+        "Output of H(0) and O(0) can only be used with redox.");
  }
+
+  this->with_h0_o0 = with_h0_o0;
+  this->with_redox = with_redox;
+
+  PhreeqcMatrix pqc_mat =
+      PhreeqcMatrix(database, script, with_h0_o0, with_redox);
+
+  this->transport_names = pqc_mat.getMatrixTransported();
+
+  Rcpp::Function unique_R("unique");
+  Rcpp::Function sort_R("sort");
+
+  std::vector<int> unique_ids = Rcpp::as<std::vector<int>>(
+      unique_R(Rcpp::IntegerVector(grid_def.begin(), grid_def.end())));
+
+  this->initial_grid = expandGrid(pqc_mat, unique_ids, grid_def);
+
+  const auto pqc_raw_dumps = replaceRawKeywordIDs(pqc_mat.getDumpStringsPQI());
+
+  for (const auto &id : unique_ids) {
+    this->pqc_ids.push_back(id);
+  }
+
+  return pqc_mat;
+  // this->phreeqc_mat = pqcScriptToGrid(phreeqc, R);
+  // this->initial_grid = matToGrid(R, this->phreeqc_mat, grid_def);
+
+  // const uint32_t solution_count = getSolutionCount(phreeqc,
+  // this->initial_grid);
+
+  // std::vector<std::string> colnames =
+  //     Rcpp::as<std::vector<std::string>>(this->initial_grid.names());
+
+  // this->transport_names = std::vector<std::string>(
+  //     colnames.begin() + 1,
+  //     colnames.begin() + 1 + solution_count); // skip ID
+
+  // std::map<int, std::string> pqc_raw_dumps;
+
+  // pqc_raw_dumps = replaceRawKeywordIDs(phreeqc->raw_dumps());
+
+  // this->pqc_ids =
+  //     Rcpp::as<std::vector<int>>(unique_R(this->initial_grid["ID"]));
+
+  // for (const auto &id : this->pqc_ids) {
+  //   this->pqc_scripts.push_back(pqc_raw_dumps[id]);
+  //   // this->pqc_exchanger.push_back(phreeqc->getExchanger(id));
+  //   // this->pqc_kinetics.push_back(phreeqc->getKineticsNames(id));
+  //   // this->pqc_equilibrium.push_back(phreeqc->getEquilibriumNames(id));
+  //   // this->pqc_surface_comps.push_back(phreeqc->getSurfaceCompNames(id));
+  //   //
+  //   this->pqc_surface_charges.push_back(phreeqc->getSurfaceChargeNames(id));
+  // }
 }

-} // namespace poet
+} // namespace poet
--- a/src/Init/InitialList.cpp
+++ b/src/Init/InitialList.cpp
@ -1,5 +1,6 @@
 #include "InitialList.hpp"
 #include "DataStructures/NamedVector.hpp"
+#include "PhreeqcMatrix.hpp"
 #include <Rcpp/internal/wrap.h>
 #include <Rcpp/iostream/Rstreambuf.h>
 #include <Rcpp/proxy/ProtectedProxy.h>
@ -10,8 +11,8 @@

 namespace poet {
 void InitialList::initializeFromList(const Rcpp::List &setup) {
-  prepareGrid(setup[grid_key]);
-  initDiffusion(setup[diffusion_key]);
+  PhreeqcMatrix phreeqc = prepareGrid(setup[grid_key]);
+  initDiffusion(setup[diffusion_key], phreeqc);
  initChemistry(setup[chemistry_key]);
 }

@ -53,26 +54,30 @@ void InitialList::importList(const Rcpp::List &setup, bool minimal) {
  }
  this->database =
      Rcpp::as<std::string>(setup[static_cast<int>(ExportList::CHEM_DATABASE)]);
+  this->pqc_script = Rcpp::as<std::string>(
+      setup[static_cast<int>(ExportList::CHEM_PQC_SCRIPT)]);
+  this->with_h0_o0 =
+      Rcpp::as<bool>(setup[static_cast<int>(ExportList::CHEM_PQC_WITH_H0_O0)]);
+  this->with_redox =
+      Rcpp::as<bool>(setup[static_cast<int>(ExportList::CHEM_PQC_WITH_REDOX)]);
  this->field_header = Rcpp::as<std::vector<std::string>>(
      setup[static_cast<int>(ExportList::CHEM_FIELD_HEADER)]);
-  this->pqc_scripts = Rcpp::as<std::vector<std::string>>(
-      setup[static_cast<int>(ExportList::CHEM_PQC_SCRIPTS)]);
  this->pqc_ids = Rcpp::as<std::vector<int>>(
      setup[static_cast<int>(ExportList::CHEM_PQC_IDS)]);
-  this->pqc_solutions = Rcpp::as<std::vector<std::string>>(
-      setup[static_cast<int>(ExportList::CHEM_PQC_SOLUTIONS)]);
-  this->pqc_solution_primaries = Rcpp::as<std::vector<std::string>>(
-      setup[static_cast<int>(ExportList::CHEM_PQC_SOLUTION_PRIMARY)]);
-  this->pqc_exchanger =
-      Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_EXCHANGER)]);
-  this->pqc_kinetics =
-      Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_KINETICS)]);
-  this->pqc_equilibrium =
-      Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_EQUILIBRIUM)]);
-  this->pqc_surface_comps =
-      Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_SURFACE_COMPS)]);
-  this->pqc_surface_charges =
-      Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_SURFACE_CHARGES)]);
+  //   this->pqc_solutions = Rcpp::as<std::vector<std::string>>(
+  //       setup[static_cast<int>(ExportList::CHEM_PQC_SOLUTIONS)]);
+  //   this->pqc_solution_primaries = Rcpp::as<std::vector<std::string>>(
+  //       setup[static_cast<int>(ExportList::CHEM_PQC_SOLUTION_PRIMARY)]);
+  //   this->pqc_exchanger =
+  //       Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_EXCHANGER)]);
+  //   this->pqc_kinetics =
+  //       Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_KINETICS)]);
+  //   this->pqc_equilibrium =
+  //       Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_EQUILIBRIUM)]);
+  //   this->pqc_surface_comps =
+  //       Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_SURFACE_COMPS)]);
+  //   this->pqc_surface_charges =
+  //       Rcpp::List(setup[static_cast<int>(ExportList::CHEM_PQC_SURFACE_CHARGES)]);

  this->dht_species = NamedVector<uint32_t>(
      setup[static_cast<int>(ExportList::CHEM_DHT_SPECIES)]);
@ -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)] =
-      Rcpp::wrap(this->pqc_solutions);
-  out[static_cast<int>(ExportList::CHEM_PQC_SOLUTION_PRIMARY)] =
-      Rcpp::wrap(this->pqc_solution_primaries);
-  out[static_cast<int>(ExportList::CHEM_PQC_EXCHANGER)] =
-      Rcpp::wrap(this->pqc_exchanger);
-  out[static_cast<int>(ExportList::CHEM_PQC_KINETICS)] =
-      Rcpp::wrap(this->pqc_kinetics);
-  out[static_cast<int>(ExportList::CHEM_PQC_EQUILIBRIUM)] =
-      Rcpp::wrap(this->pqc_equilibrium);
-  out[static_cast<int>(ExportList::CHEM_PQC_SURFACE_COMPS)] =
-      Rcpp::wrap(this->pqc_surface_comps);
-  out[static_cast<int>(ExportList::CHEM_PQC_SURFACE_CHARGES)] =
-      Rcpp::wrap(this->pqc_surface_charges);
+  //   out[static_cast<int>(ExportList::CHEM_PQC_SOLUTIONS)] =
+  //       Rcpp::wrap(this->pqc_solutions);
+  //   out[static_cast<int>(ExportList::CHEM_PQC_SOLUTION_PRIMARY)] =
+  //       Rcpp::wrap(this->pqc_solution_primaries);
+  //   out[static_cast<int>(ExportList::CHEM_PQC_EXCHANGER)] =
+  //       Rcpp::wrap(this->pqc_exchanger);
+  //   out[static_cast<int>(ExportList::CHEM_PQC_KINETICS)] =
+  //       Rcpp::wrap(this->pqc_kinetics);
+  //   out[static_cast<int>(ExportList::CHEM_PQC_EQUILIBRIUM)] =
+  //       Rcpp::wrap(this->pqc_equilibrium);
+  //   out[static_cast<int>(ExportList::CHEM_PQC_SURFACE_COMPS)] =
+  //       Rcpp::wrap(this->pqc_surface_comps);
+  //   out[static_cast<int>(ExportList::CHEM_PQC_SURFACE_CHARGES)] =
+  //       Rcpp::wrap(this->pqc_surface_charges);
  out[static_cast<int>(ExportList::CHEM_DHT_SPECIES)] = this->dht_species;
  out[static_cast<int>(ExportList::CHEM_INTERP_SPECIES)] =
      Rcpp::wrap(this->interp_species);
--- a/src/Init/InitialList.hpp
+++ b/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>
@ -35,7 +33,7 @@ public:
  void importList(const Rcpp::List &setup, bool minimal = false);
  Rcpp::List exportList();

-  Field getInitialGrid() const { return Field(this->initial_grid); }  
+  Field getInitialGrid() const { return Field(this->initial_grid); }

 private:
  RInside &R;
@ -53,22 +51,17 @@ private:
    DIFFU_ALPHA_X,
    DIFFU_ALPHA_Y,
    CHEM_DATABASE,
-    CHEM_FIELD_HEADER,
-    CHEM_PQC_SCRIPTS,
+    CHEM_PQC_SCRIPT,
+    CHEM_PQC_WITH_H0_O0,
+    CHEM_PQC_WITH_REDOX,
    CHEM_PQC_IDS,
-    CHEM_PQC_SOLUTIONS,
-    CHEM_PQC_SOLUTION_PRIMARY,
-    CHEM_PQC_EXCHANGER,
-    CHEM_PQC_KINETICS,
-    CHEM_PQC_EQUILIBRIUM,
-    CHEM_PQC_SURFACE_COMPS,
-    CHEM_PQC_SURFACE_CHARGES,
+    CHEM_FIELD_HEADER,
    CHEM_DHT_SPECIES,
    CHEM_INTERP_SPECIES,
    CHEM_HOOKS,
    AI_SURROGATE_INPUT_SCRIPT,
    ENUM_SIZE // Hack: Last element of the enum to show enum size
-    };
+  };

  // Grid members
  static constexpr const char *grid_key = "Grid";
@ -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",
-                           "pqc_db_file",    "grid_def",    "grid_size",
-                           "constant_cells", "porosity"};
+      GridMembersString = {"pqc_in_string", "pqc_in_file",   "pqc_with_redox",
+                           "pqc_wth_h0_o0", "pqc_db_string", "pqc_db_file",
+                           "grid_def",      "grid_size",     "constant_cells",
+                           "porosity"};

  constexpr const char *GRID_MEMBER_STR(GridMembers member) const {
    return GridMembersString[static_cast<std::size_t>(member)];
  }

-  std::unique_ptr<PhreeqcInit> phreeqc;
+  // std::unique_ptr<PhreeqcMatrix> pqc_mat;

-  void prepareGrid(const Rcpp::List &grid_input);
+  PhreeqcMatrix prepareGrid(const Rcpp::List &grid_input);

  std::uint8_t dim{0};

@ -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,21 +185,16 @@ private:
  std::vector<std::string> field_header;

  std::string database;
-  std::vector<std::string> pqc_scripts;
+  std::string pqc_script;
+  bool with_h0_o0{false};
+  bool with_redox{false};
+  // std::vector<std::string> pqc_scripts;
  std::vector<int> pqc_ids;

-  std::vector<std::string> pqc_solutions;
-  std::vector<std::string> pqc_solution_primaries;
-  Rcpp::List pqc_exchanger;
-  Rcpp::List pqc_kinetics;
-  Rcpp::List pqc_equilibrium;
-  Rcpp::List pqc_surface_comps;
-  Rcpp::List pqc_surface_charges;
-
  NamedVector<std::uint32_t> dht_species;

  NamedVector<std::uint32_t> interp_species;
-  
+
  // Path to R script that the user defines in the input file
  std::string ai_surrogate_input_script;

@ -227,10 +217,12 @@ public:
    // std::vector<std::string> field_header;

    std::string database;
-    // std::vector<std::string> pqc_scripts;
-    // std::vector<int> pqc_ids;
+    std::string pqc_script;
+    bool with_h0_o0;
+    bool with_redox;
+    std::vector<int> pqc_ids;

-    std::map<int, POETConfig> pqc_config;
+    // std::map<int, std::string> pqc_input;

    // std::vector<std::string> pqc_sol_order;

--- a/src/Transport/DiffusionModule.cpp
+++ b/src/Transport/DiffusionModule.cpp
@ -1,97 +1,48 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
 /*
-** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
-** Potsdam)
-**
-** Copyright (C) 2018-2022 Marco De Lucia, Max Luebke (GFZ Potsdam)
-**
-** Copyright (C) 2023-2024 Marco De Lucia (GFZ Potsdam), Max Luebke (University
-** of Potsdam)
-**
-** POET is free software; you can redistribute it and/or modify it under the
-** terms of the GNU General Public License as published by the Free Software
-** Foundation; either version 2 of the License, or (at your option) any later
-** version.
-**
-** POET is distributed in the hope that it will be useful, but WITHOUT ANY
-** WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
-** A PARTICULAR PURPOSE. See the GNU General Public License for more details.
-**
-** You should have received a copy of the GNU General Public License along with
-** this program; if not, write to the Free Software Foundation, Inc., 51
-** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
-*/
+ * DiffusionModule.cpp - Coupling module betweeen POET and tug
+ * Copyright (C) 2018-2025 Max Luebke (University of Potsdam)
+ */

 #include "DiffusionModule.hpp"

 #include "Base/Macros.hpp"
 #include "Init/InitialList.hpp"

-#include <Eigen/src/Core/Map.h>
-#include <Eigen/src/Core/Matrix.h>
-#include <Eigen/src/Core/util/Constants.h>
 #include <chrono>
-
-#include "tug/Boundary.hpp"
-#include "tug/Grid.hpp"
-#include "tug/Simulation.hpp"
-
 #include <string>
-#include <vector>
+#include <tug/Boundary.hpp>
+#include <tug/Diffusion.hpp>

 using namespace poet;

-static inline std::vector<TugType>
-MatrixToVec(const Eigen::MatrixX<TugType> &mat) {
-  std::vector<TugType> vec(mat.rows() * mat.cols());
-
-  for (std::uint32_t i = 0; i < mat.cols(); i++) {
-    for (std::uint32_t j = 0; j < mat.rows(); j++) {
-      vec[j * mat.cols() + i] = mat(j, i);
-    }
-  }
-
-  return vec;
-}
-
-static inline Eigen::MatrixX<TugType>
-VecToMatrix(const std::vector<TugType> &vec, std::uint32_t n_rows,
-            std::uint32_t n_cols) {
-  Eigen::MatrixX<TugType> mat(n_rows, n_cols);
-
-  for (std::uint32_t i = 0; i < n_cols; i++) {
-    for (std::uint32_t j = 0; j < n_rows; j++) {
-      mat(j, i) = vec[j * n_cols + i];
-    }
-  }
-
-  return mat;
-}
-
-// static constexpr double ZERO_MULTIPLICATOR = 10E-14;
-
 void DiffusionModule::simulate(double requested_dt) {
-  MSG("Starting diffusion ...");
+  // MSG("Starting diffusion ...");
  const auto start_diffusion_t = std::chrono::high_resolution_clock::now();

  const auto &n_rows = this->param_list.n_rows;
  const auto &n_cols = this->param_list.n_cols;

-  tug::Grid<TugType> grid(param_list.n_rows, param_list.n_cols);
-  tug::Boundary boundary(grid);
-
-  grid.setDomain(param_list.s_rows, param_list.s_cols);
-
-  tug::Simulation<TugType> sim(grid, boundary);
-  sim.setIterations(1);
-
  for (const auto &sol_name : this->param_list.transport_names) {
-    auto &species_conc = this->transport_field[sol_name];
+#if defined(POET_TUG_BTCS)
+    tug::Diffusion<TugType> diffusion_solver(
+        this->transport_field[sol_name].data(), n_rows, n_cols);
+#elif defined(POET_TUG_FTCS)
+    tug::Diffusion<TugType, tug::FTCS_APPROACH> diffusion_solver(
+        this->transport_field[sol_name].data(), n_rows, n_cols);
+#else
+#error "No valid approach defined"
+#endif

-    Eigen::MatrixX<TugType> conc = VecToMatrix(species_conc, n_rows, n_cols);
-    Eigen::MatrixX<TugType> alpha_x =
-        VecToMatrix(this->param_list.alpha_x[sol_name], n_rows, n_cols);
-    Eigen::MatrixX<TugType> alpha_y =
-        VecToMatrix(this->param_list.alpha_y[sol_name], n_rows, n_cols);
+    tug::RowMajMatMap<TugType> alpha_x_map(
+        this->param_list.alpha_x[sol_name].data(), n_rows, n_cols);
+    tug::RowMajMatMap<TugType> alpha_y_map(
+        this->param_list.alpha_y[sol_name].data(), n_rows, n_cols);
+
+    diffusion_solver.setAlphaX(alpha_x_map);
+    diffusion_solver.setAlphaY(alpha_y_map);
+
+    auto &boundary = diffusion_solver.getBoundaryConditions();

    boundary.deserialize(this->param_list.boundaries[sol_name]);

@ -99,14 +50,8 @@ void DiffusionModule::simulate(double requested_dt) {
      boundary.setInnerBoundaries(this->param_list.inner_boundaries[sol_name]);
    }

-    grid.setAlpha(alpha_x, alpha_y);
-    grid.setConcentrations(conc);
-
-    sim.setTimestep(requested_dt);
-
-    sim.run();
-
-    species_conc = MatrixToVec(grid.getConcentrations());
+    diffusion_solver.setTimestep(requested_dt);
+    diffusion_solver.run();
  }

  const auto end_diffusion_t = std::chrono::high_resolution_clock::now();
--- a/src/initializer.cpp
+++ b/src/initializer.cpp
@ -34,9 +34,11 @@ int main(int argc, char **argv) {
               "input script")
      ->default_val(false);

-  bool asRDS;
-  app.add_flag("-r, --rds", asRDS, "Save output as .rds file instead of .qs")
-      ->default_val(false);
+  bool asRDS{false};
+  app.add_flag("-r, --rds", asRDS, "Save output as .rds")->default_val(false);
+
+  bool asQS{false};
+  app.add_flag("-q, --qs", asQS, "Save output as .qs")->default_val(false);

  CLI11_PARSE(app, argc, argv);

@ -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) {
--- a/src/poet.cpp
+++ b/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;
-    MSG("Next time step is " + std::to_string(dt) + " [s]");
+    std::cout << std::endl;

    /* displaying iteration number, with C++ and R iterator */
-    MSG("Going through iteration " + std::to_string(iter));
+    MSG("Going through iteration " + std::to_string(iter) + "/" +
+        std::to_string(maxiter));
+
+    MSG("Current time step is " + std::to_string(dt));

    /* run transport */
    diffusion.simulate(dt);

    chem.getField().update(diffusion.getField());

-    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)");

@ -333,8 +345,8 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, const RuntimeParameters &params,
      MSG("AI TempField");
      std::vector<std::vector<double>> RTempField =
          R.parseEval("set_valid_predictions(predictors,\
-                                                                                       aipreds,\
-                                                                                       validity_vector)");
+                       aipreds,\
+                       validity_vector)");

      MSG("AI Set Field");
      Field predictions_field =
@ -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);
--- a/src/poet.hpp.in
+++ b/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;
 };
Author	SHA1	Message	Date
Marco De Lucia	58175dd8e2	Merge branch 'mdl/defaultstorelease' into 'main' DEFAULT_BUILD_TYPE set to Release See merge request naaice/poet!60	2025-12-06 22:37:17 +01:00
Marco De Lucia	c8a779ca58	set DEFAULT_BUILD_TYPE to "Release"; defined properties so one cycles through when using ccmake	2025-11-14 19:19:07 +01:00
Marco De Lucia	14c08460ad	Merge branch 'mdl/readme' into 'main' Mdl/readme Closes #22 See merge request naaice/poet!59	2025-11-14 15:31:46 +01:00
Marco De Lucia	0dc9352459	Updated README - domain change, litephreeqc	2025-11-14 15:28:25 +01:00
Max Lübke	9cd487cc15	Merge branch 'mdl/litephreeqc' into 'main' Chore: New URL and name for litephreeqc in .gitmodules See merge request naaice/poet!58	2025-11-13 16:45:53 +01:00
Max Lübke	04ac4b8b7b	ci(gitlab): allow mirror push job to fail	2025-11-13 16:39:06 +01:00
Max Lübke	2c813b9e24	build: set path for submodule	2025-11-13 16:35:57 +01:00
Max Lübke	1b55f5e651	chore: use git mv to set new path for litephreeqc	2025-11-13 16:31:47 +01:00
Max Lübke	746ae15645	chore: rollback submodule path	2025-11-13 16:30:55 +01:00
Marco De Lucia	9066ae942b	Chore: New URL and name for litephreeqc in .gitmodules	2025-11-13 16:22:21 +01:00
Marco De Lucia	22f5f6f21d	chore: add CITATIONS.cff	2025-10-24 11:29:57 +02:00
Marco De Lucia	cbe0712766	Merge branch '19-cmake-requirements-missmatch-with-tug' into 'main' Resolve "CMake requirements missmatch with TUG." Closes #19 See merge request naaice/poet!57	2025-10-24 10:33:41 +02:00
Max Lübke	21c606f5b5	docs(license): add SPDX identifier to source files	2025-10-24 09:41:25 +02:00
Max Lübke	c8bedaac28	refactor(diffusion): integrate tug::Diffusion solver	2025-10-24 09:41:03 +02:00
Max Lübke	5f2b0cc942	build(deps): update cli11 to v2.5.0	2025-10-24 09:38:26 +02:00
Max Lübke	0d307c790f	build(deps): update minimum CMake to 3.20 and tug submodule	2025-10-24 09:38:05 +02:00
Max Lübke	fb1f053055	Merge branch 'ml/mirror-to-github' into 'main' ci(github): add job to push to GitHub remote See merge request naaice/poet!56	2025-10-15 09:58:35 +02:00
Max Lübke	86a683f07c	ci(github): add job to push to GitHub remote	2025-10-15 09:57:05 +02:00
Max Lübke	490a28903b	Merge branch 'chore' into 'main' chore: update iphreeqc submodule See merge request naaice/poet!54	2025-09-23 12:21:14 +02:00
Max Lübke	0317ddbf8f	chore: update iphreeqc submodule	2025-09-23 12:19:45 +02:00
Max Lübke	92d6044e97	Merge branch 'ml/add-build-type-to-readme' into 'main' Add CMake Build Type to Readme See merge request naaice/poet!53	2025-05-20 14:29:00 +02:00
Max Lübke	77ccd11220	Add CMake Build Type to Readme	2025-05-20 14:28:45 +02:00
Marco De Lucia	ccd54aaa09	Merge branch '17-new-barite-benchmarks-for-fgcs-journal' into 'main' Resolve "New barite benchmarks for FGCS journal" Closes #17 See merge request naaice/poet!44	2025-05-19 12:50:30 +02:00
Max Lübke	66d908b548	fix: Correct boolean conversion in DHT hook call	2025-03-12 10:35:23 +01:00
Hannes Martin Signer	a4bf2a13f9	Merge branch 'ml/add-tug-solver-option' into 'main' Add option to toggle tug numerical solver approach See merge request naaice/poet!50	2025-01-29 16:17:52 +01:00
Max Luebke	85e93955bc	Add option to toggle tug numerical solver approach	2025-01-29 15:53:32 +01:00
Max Lübke	5a7779e8de	refactor: clean up CMakeLists.txt by removing unnecessary target_sources indentation	2025-01-09 12:18:13 +01:00
Max Lübke	bcab85c331	feat: add index tracking and flushing functionality to DHT	2025-01-09 12:18:05 +01:00
Max Lübke	dffbec674c	Merge branch '17-new-barite-benchmarks-for-fgcs-journal' of git.gfz-potsdam.de:naaice/poet into 17-new-barite-benchmarks-for-fgcs-journal	2024-12-20 10:12:51 +01:00
Max Lübke	2df3f9f34a	[wip] add input validation for work_package in DHT_Wrapper and InterpolationModule	2024-12-20 10:12:19 +01:00
Max Lübke	cefea926cf	fix: update grid dimensions and refine species definitions in barite_fgcs_2.R	2024-12-20 10:08:07 +01:00
Marco De Lucia	b93cdeed83	Fancier benchmark for fgcs paper	2024-12-20 10:08:07 +01:00
Max Lübke	8b7e3e983e	Merge branch 'hotfix-qs-option' into 'main' fix: remove duplicate flag for saving output as .qs file in parseInitValues See merge request naaice/poet!49	2024-12-20 10:07:21 +01:00
Max Lübke	a14095622f	fix: remove duplicate flag for saving output as .qs file in parseInitValues	2024-12-20 10:05:49 +01:00
Marco De Lucia	554d501ed0	Fancier benchmark for fgcs paper	2024-12-20 09:24:19 +01:00
Max Lübke	ab01dbc0a7	Merge branch 'ml/fix-interpolation' into 'main' Fix Interpolation and add optional features to POET benchmark description See merge request naaice/poet!48	2024-12-20 09:23:44 +01:00
Max Lübke	eb384b48c0	Merge branch 'main' into 'ml/fix-interpolation' # Conflicts: # src/initializer.cpp # src/poet.cpp # src/poet.hpp.in	2024-12-20 09:22:55 +01:00
Max Lübke	e34c2ef464	fix: initialize default values for RuntimeParameters in poet.hpp.in	2024-12-20 09:21:18 +01:00
Max Lübke	d5a70a20bb	fix: update dht_species values for barite and celestite in benchmark scripts	2024-12-20 09:21:18 +01:00
Max Lübke	a718c3efaf	Revert "feat: make output of H(0) and O(0) optional" This reverts commit bdb5ce944fd001cb36c8369d8c792f2f4a54c7be.	2024-12-20 09:21:18 +01:00
Max Lübke	5664153cf9	fix: handle zero distance case in inverseDistanceWeighting to prevent division by zero	2024-12-20 09:21:18 +01:00
Max Lübke	6ab404eefd	feat: add with_h0_o0 and with_redox parameters to ChemistryModule and related classes	2024-12-20 09:21:18 +01:00
Max Luebke	e5588a2ee9	feat: add support for redox in PhreeqcMatrix initialization	2024-12-20 09:21:18 +01:00
Max Lübke	6fb226887e	feat: enhance LookupKey with const isnan method and comparison operator	2024-12-20 09:21:18 +01:00
Max Lübke	521937a527	refactor: streamline WorkPackage constructor and improve DHT_Wrapper::fillDHT logic	2024-12-20 09:21:17 +01:00
Max Luebke	d661da13bd	feat: add dht_snaps and dht_out_dir parameters to ChemistryModule and main function	2024-12-20 09:21:17 +01:00
Max Luebke	06728bc28b	fix: initialize includeH0O0 to false in command line options	2024-12-20 09:21:17 +01:00
Max Luebke	c5b0ae201c	feat: make output of H(0) and O(0) optional	2024-12-20 09:21:17 +01:00
Max Lübke	ac96f24a33	feat: add has_het_ids parameter to DHT initialization and related functions	2024-12-20 09:21:17 +01:00
Max Lübke	9b9fd898b7	fix: use const reference for to_calc in InterpolationModule to improve performance	2024-12-20 09:21:17 +01:00
Max Lübke	536ebdd351	fix: remove double initialization of dht_enabled in ChemistryModule	2024-12-20 09:21:17 +01:00
Max Lübke	f6b4ce017a	feat: implement caching for interpolation calculations and add NaN handling	2024-12-20 09:21:16 +01:00
Max Luebke	43d2a846c7	[wip] fix: add base_totals to SurrogateSetup	2024-12-20 09:21:16 +01:00
Marco De Lucia	dd9cc5e59f	reverting <format> since gcc < 13 does not support it	2024-12-20 09:21:16 +01:00
Marco De Lucia	2e3de84a90	Update references to .qs2 in README	2024-12-20 09:21:16 +01:00
Marco De Lucia	3b47b8ac1f	Update README: qs2 as default output format, gfz.de everywhere	2024-12-20 09:21:16 +01:00
Marco De Lucia	0437670e05	Less and more informative stdout messages	2024-12-20 09:21:16 +01:00
Marco De Lucia	9f89edd492	Added qs2 as new default format	2024-12-20 09:21:15 +01:00
Max Lübke	347b60c3f3	Merge branch 'ml/gitlab-ci' into 'main' Refactor CI pipeline by removing the build stage and adjusting test job dependencies See merge request naaice/poet!47	2024-12-20 09:06:18 +01:00
Max Lübke	c36fe346b3	Refactor CI pipeline by removing the build stage and adjusting test job dependencies	2024-12-20 09:04:30 +01:00
Max Lübke	f4f4ea2be2	Merge branch 'ml/fix-license' into 'main' Revert "No code changes made." See merge request naaice/poet!46	2024-12-20 09:00:27 +01:00
Max Lübke	60c4a15c9f	Update LICENSE to include updated copyright notice and additional guidelines for program distribution	2024-12-20 08:58:40 +01:00
Max Lübke	0c22f17df7	Revert "No code changes made." This reverts commit f60b7cf7fcb33e9ff55e1c8050775217d9d37f1a.	2024-12-20 08:57:51 +01:00
Max Lübke	824442b315	Merge branch 'ml/fix-license' into 'main' Fix license parsing error of GitLab? See merge request naaice/poet!45	2024-12-20 08:55:43 +01:00
Max Lübke	f60b7cf7fc	No code changes made.	2024-12-20 08:46:51 +01:00
Max Lübke	d121585894	Merge branch 'mdl/fgcsbench' into 'main' Mdl/fgcsbench See merge request naaice/poet!43	2024-12-13 12:19:57 +01:00
Marco De Lucia	03385d3caf	reverting <format> since gcc < 13 does not support it	2024-12-13 12:14:55 +01:00
Marco De Lucia	2478518d73	Update references to .qs2 in README	2024-12-13 12:14:55 +01:00
Marco De Lucia	9a787b793b	Update README: qs2 as default output format, gfz.de everywhere	2024-12-13 12:14:55 +01:00
Marco De Lucia	a519edd102	Less and more informative stdout messages	2024-12-13 12:14:55 +01:00
Marco De Lucia	61d3a8a227	Added qs2 as new default format	2024-12-13 12:14:55 +01:00
Hannes Martin Signer	b29be5fd6d	Merge branch 'update-scheme' into 'main' update poet scheme See merge request naaice/poet!42	2024-12-10 15:03:54 +01:00
Hannes Martin Signer	b19cdef894	Update 2 files - /docs/POET_scheme.svg - /docs/POET.drawio	2024-12-10 15:02:45 +01:00
Max Lübke	b67c1fd36d	Merge branch 'update-poet-scheme' into 'main' update POET scheme See merge request naaice/poet!41	2024-12-10 11:15:48 +01:00
Hannes Martin Signer	2dba14d3bd	update POET scheme	2024-12-10 11:14:01 +01:00
Max Lübke	5c92a2156f	Merge branch 'ml/iphreeqc-v0.3' into 'main' Update IPhreeqc Submodule to v0.3 See merge request naaice/poet!40	2024-12-04 11:35:27 +01:00
Max Lübke	6c660557fd	refactor: Replace PhreeqcEngine instances with PhreeqcRunner for improved chemistry processing	2024-12-04 10:34:13 +00:00
Max Lübke	db7a2ad2ce	build: Enhance Dockerfile for improved environment setup and dependency management	2024-12-04 08:44:55 +00:00
Max Lübke	0d9162fb66	Merge branch 'ml/misc' into 'main' Various changes to build files (CMake & devcontainer) See merge request naaice/poet!38	2024-11-07 20:57:45 +01:00
Max Lübke	716454e115	build: Update Dockerfile and devcontainer.json for enhanced environment setup	2024-11-07 19:56:04 +00:00
Max Lübke	287eda880a	build: Refactor R runtime detection in CMake to ensure required parameters are set	2024-11-07 19:56:04 +00:00
Max Lübke	66098aab40	Merge branch 'ml/new-iphreeqc-api' into 'main' Update to IPhreeqc v3.8.2 See merge request naaice/poet!37	2024-11-07 14:23:58 +01:00
Max Lübke	cdbf344329	fix: Update getSolutionNames call to remove unnecessary argument	2024-11-07 14:22:58 +01:00
Max Lübke	45ea77ae0f	Update subproject commit in ext/iphreeqc	2024-11-07 14:22:51 +01:00
				`@ -0,0 +1,2 @@`
				`* Refer to the LaTeX file (and pdf) for more information`
				`@ -1 +0,0 @@`
				`Subproject commit e6e5e0d5156c093241a53e6ce074ef346d64ae26`
				`@ -0,0 +1 @@`
				`Subproject commit 953c752431d2b2758268083f407f943843efc7ad`