max/poet
1
0
Fork 0
mirror of https://git.gfz-potsdam.de/naaice/poet.git synced 2025-12-16 04:48:23 +01:00
Code Issues Packages Projects Releases Wiki Activity

Compare commits

base: max:main
Branches Tags
max:main
max:naa-naaice
max:naaice-baseline
max:ml/add_real_barite_benchmark
max:control-loop-proto1
max:control-loop-proto2
max:readme
max:wip-selected-output
max:refac/hd5
max:tmp-fgcs
max:serialization-change
max:naaice
max:mdl/fmt
max:ml/field_memory
max:0.4
max:archive/iccs2025
max:archive/hans_ai_surrogate
max:archive/hans_k_means
max:v0.3
max:v0.2.1
max:archive/ns-daos-diffusion-2
max:archive/ns-daos-diffusion-1
max:v0.2
max:v0.1
..
compare: max:v0.3
Branches Tags
max:naa-naaice
max:naaice-baseline
max:ml/add_real_barite_benchmark
max:control-loop-proto1
max:control-loop-proto2
max:main
max:readme
max:wip-selected-output
max:refac/hd5
max:tmp-fgcs
max:serialization-change
max:naaice
max:mdl/fmt
max:ml/field_memory
max:0.4
max:archive/iccs2025
max:archive/hans_ai_surrogate
max:archive/hans_k_means
max:v0.3
max:v0.2.1
max:archive/ns-daos-diffusion-2
max:archive/ns-daos-diffusion-1
max:v0.2
max:v0.1

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

69 changed files with 2148 additions and 3354 deletions
Show Stats Expand all files Collapse all files

122
.devcontainer/Dockerfile
View File

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

3
.devcontainer/devcontainer.json
View File

@ -21,8 +21,7 @@
},
// in case you want to push/pull from remote repositories using ssh
"mounts": [
"source=${localEnv:HOME}/.ssh,target=/home/vscode/.ssh,type=bind,consistency=cached",
"source=${localEnv:HOME}/.gitconfig,target=/home/vscode/.gitconfig,type=bind,consistency=cached"
"source=${localEnv:HOME}/.ssh,target=/home/vscode/.ssh,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"

38
.gitlab-ci.yml
View File

@ -20,6 +20,7 @@ image: git.gfz-potsdam.de:5000/naaice/poet:ci
stages: # List of stages for jobs, and their order of execution
- release
- build
- test
variables:
@ -27,11 +28,26 @@ variables:
SOURCE_ARCHIVE_NAME: 'poet_${CI_COMMIT_TAG}_sources.tar.gz'
CHANGELOG_FILE: 'commit_changelog.md'
test: # This job runs in the build stage, which runs first.
stage: test
build-poet: # This job runs in the build stage, which runs first.
stage: build
script:
- mkdir -p build && cd build
- cmake -DPOET_ENABLE_TESTING=ON -DPOET_PREPROCESS_BENCHS=OFF -DCMAKE_BUILD_TYPE=Release ..
- 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
- make -j$(nproc) check
pages:
@ -49,22 +65,6 @@ 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

6
.gitmodules vendored
View File

@ -2,6 +2,6 @@
path = ext/tug
url = ../tug.git
[submodule "ext/litephreeqc"]
path = ext/litephreeqc
url = ../litephreeqc.git
[submodule "ext/iphreeqc"]
path = ext/iphreeqc
url = ../iphreeqc.git

51
CITATION.cff
View File

@ -1,51 +0,0 @@
# This CITATION.cff file was generated with cffinit.
# Visit https://bit.ly/cffinit to generate yours today!
cff-version: 1.2.0
title: 'POET: POtsdam rEactive Transport'
message: >-
If you use this software, please cite it using these
metadata.
type: software
authors:
- given-names: Max
family-names: Lübke
email: mluebke@uni-potsdam.de
affiliation: University of Potsdam
orcid: 'https://orcid.org/0009-0008-9773-3038'
- given-names: Marco
family-names: De Lucia
email: delucia@gfz.de
affiliation: GFZ Helmholtz Centre for Geosciences
orcid: 'https://orcid.org/0000-0002-1186-4491'
- given-names: Alexander
family-names: Lindemann
- given-names: Hannes
family-names: Signer
email: signer@uni-potsdam.de
orcid: 'https://orcid.org/0009-0000-3058-8472'
- given-names: Bettina
family-names: Schnor
email: schnor@cs.uni-potsdam.de
affiliation: University of Potsdam
orcid: 'https://orcid.org/0000-0001-7369-8057'
- given-names: Hans
family-names: Straile
identifiers:
- type: doi
value: 10.5194/gmd-14-7391-2021
description: >-
POET (v0.1): speedup of many-core parallel reactive
transport simulations with fast DHT lookups
repository-code: 'https://git.gfz-potsdam.de/naaice/poet'
abstract: >-
Massively parallel reactive transport simulator exploring
acceleration strategies such as embedding of AI/ML and
cache of results in Distributed Hash Tables. Developed in
cooperation with computer scientists of University of
Potsdam.
keywords:
- Reactive Transport
- Geochemistry
- AI/ML Surrogate Modelling
license: GPL-2.0

52
CMake/FindRRuntime.cmake
View File

@ -1,27 +1,30 @@
# prepare R environment (Rcpp + RInside)
find_program(R_EXE "R"
REQUIRED
)
find_program(R_EXE "R")
# search for R executable, R header file and library path
execute_process(
COMMAND ${R_EXE} RHOME
OUTPUT_VARIABLE R_ROOT_DIR
OUTPUT_STRIP_TRAILING_WHITESPACE
)
if(R_EXE)
execute_process(
COMMAND ${R_EXE} RHOME
OUTPUT_VARIABLE R_ROOT_DIR
OUTPUT_STRIP_TRAILING_WHITESPACE
)
find_path(
R_INCLUDE_DIR R.h
HINTS /usr/include /usr/local/include /usr/share ${R_ROOT_DIR}/include
PATH_SUFFIXES R/include R
REQUIRED
)
find_path(
R_INCLUDE_DIR R.h
HINTS ${R_ROOT_DIR}
PATHS /usr/include /usr/local/include /usr/share
PATH_SUFFIXES include/R R/include
)
find_library(
R_LIBRARY libR.so
HINTS ${R_ROOT_DIR}/lib
REQUIRED
)
find_library(
R_LIBRARY libR.so
HINTS ${R_ROOT_DIR}/lib
)
else()
message(FATAL_ERROR "No R runtime found!")
endif()
mark_as_advanced(R_INCLUDE_DIR R_LIBRARY R_EXE)
set(R_LIBRARIES ${R_LIBRARY})
set(R_INCLUDE_DIRS ${R_INCLUDE_DIR})
@ -42,6 +45,8 @@ find_path(R_Rcpp_INCLUDE_DIR Rcpp.h
HINTS ${RCPP_PATH}
PATH_SUFFIXES include)
mark_as_advanced(R_Rcpp_INCLUDE_DIR)
list(APPEND R_INCLUDE_DIRS ${R_Rcpp_INCLUDE_DIR})
# find RInside libraries and include path
@ -67,11 +72,16 @@ find_path(R_RInside_INCLUDE_DIR RInside.h
list(APPEND R_LIBRARIES ${R_RInside_LIBRARY})
list(APPEND R_INCLUDE_DIRS ${R_RInside_INCLUDE_DIR})
mark_as_advanced(R_RInside_LIBRARY R_RInside_INCLUDE_DIR)
# putting all together into interface library
add_library(RRuntime INTERFACE IMPORTED)
target_link_libraries(RRuntime INTERFACE ${R_LIBRARIES})
target_include_directories(RRuntime INTERFACE ${R_INCLUDE_DIRS})
set_target_properties(
RRuntime PROPERTIES
INTERFACE_LINK_LIBRARIES "${R_LIBRARIES}"
INTERFACE_INCLUDE_DIRECTORIES "${R_INCLUDE_DIRS}"
)
unset(R_LIBRARIES)
unset(R_INCLUDE_DIRS)

4
CMake/POET_Scripts.cmake
View File

@ -13,7 +13,7 @@ macro(get_POET_version)
WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}
OUTPUT_VARIABLE POET_GIT_VERSION
OUTPUT_STRIP_TRAILING_WHITESPACE)
if(NOT POET_GIT_BRANCH STREQUAL "main")
if(NOT POET_GIT_BRANCH STREQUAL "master")
set(POET_VERSION "${POET_GIT_BRANCH}/${POET_GIT_VERSION}")
else()
set(POET_VERSION "${POET_GIT_VERSION}")
@ -21,7 +21,7 @@ macro(get_POET_version)
elseif(EXISTS ${PROJECT_SOURCE_DIR}/.svn)
file(STRINGS .gitversion POET_VERSION)
else()
set(POET_VERSION "not_versioned")
set(POET_VERSION "0.1")
endif()
message(STATUS "Configuring POET version ${POET_VERSION}")

16
CMakeLists.txt
View File

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

62
LICENSE
View File

@ -1,8 +1,8 @@
GNU GENERAL PUBLIC LICENSE
Version 2, June 1991
Copyright (C) 1989, 1991 Free Software Foundation, Inc.,
<https://fsf.org/>
Copyright (C) 1989, 1991 Free Software Foundation, Inc., <http://fsf.org/>
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
Everyone is permitted to copy and distribute verbatim copies
of this license document, but changing it is not allowed.
@ -278,61 +278,3 @@ PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
POSSIBILITY OF SUCH DAMAGES.
END OF TERMS AND CONDITIONS
How to Apply These Terms to Your New Programs
If you develop a new program, and you want it to be of the greatest
possible use to the public, the best way to achieve this is to make it
free software which everyone can redistribute and change under these terms.
To do so, attach the following notices to the program. It is safest
to attach them to the start of each source file to most effectively
convey the exclusion of warranty; and each file should have at least
the "copyright" line and a pointer to where the full notice is found.
<one line to give the program's name and a brief idea of what it does.>
Copyright (C) <year> <name of author>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, see <https://www.gnu.org/licenses/>.
Also add information on how to contact you by electronic and paper mail.
If the program is interactive, make it output a short notice like this
when it starts in an interactive mode:
Gnomovision version 69, Copyright (C) year name of author
Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
This is free software, and you are welcome to redistribute it
under certain conditions; type `show c' for details.
The hypothetical commands `show w' and `show c' should show the appropriate
parts of the General Public License. Of course, the commands you use may
be called something other than `show w' and `show c'; they could even be
mouse-clicks or menu items--whatever suits your program.
You should also get your employer (if you work as a programmer) or your
school, if any, to sign a "copyright disclaimer" for the program, if
necessary. Here is a sample; alter the names:
Yoyodyne, Inc., hereby disclaims all copyright interest in the program
`Gnomovision' (which makes passes at compilers) written by James Hacker.
<signature of Moe Ghoul>, 1 April 1989
Moe Ghoul, President of Vice
This General Public License does not permit incorporating your program into
proprietary programs. If your program is a subroutine library, you may
consider it more useful to permit linking proprietary applications with the
library. If this is what you want to do, use the GNU Lesser General
Public License instead of this License.

302
README.md
View File

@ -1,108 +1,80 @@
<!--
Time-stamp: "Last modified 2023-08-02 13:55:11 mluebke"
-->
# 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](https://doi.org/10.5281/zenodo.4757913) is a coupled reactive
transport simulator implementing a parallel architecture and a fast,
original MPI-based Distributed Hash Table.
![POET's Coupling Scheme](./docs/POET_scheme.svg)
![POET's Coupling Scheme](./docs/20230720_Scheme_POET_en.svg)
## Parsed code documentiation
A parsed version of POET's documentation can be found at [Gitlab
pages](https://naaice.git-pages.gfz.de/poet).
A parsed version of POET's documentiation can be found at [Gitlab
pages](https://naaice.git-pages.gfz-potsdam.de/poet).
## External Libraries
The following external libraries are shipped with POET:
The following external header library is shipped with POET:
- **CLI11** - <https://github.com/CLIUtils/CLI11>
- **litephreeqc**: IPhreeqc
(<https://github.com/usgs-coupled/iphreeqc>) with patches from
GFZ/UP: <https://git.gfz.de/naaice/litephreeqc>
- **tug** - <https://git.gfz.de/naaice/tug>
- **argh** - https://github.com/adishavit/argh (BSD license)
- **IPhreeqc** with patches from GFZ -
https://github.com/usgs-coupled/iphreeqc -
https://git.gfz-potsdam.de/naaice/iphreeqc
- **tug** - https://git.gfz-potsdam.de/naaice/tug
## Installation
### Requirements
To compile POET you need following software to be installed:
To compile POET you need several software to be installed:
- C/C++ compiler (tested with GCC)
- MPI-Implementation (tested with OpenMPI and MVAPICH)
- CMake 3.20+
- R language and environment
- CMake 3.9+
- Eigen3 3.4+ (required by `tug`)
- *optional*: `doxygen` with `dot` bindings for documentation
- R language and environment including headers or `-dev` packages
(distro dependent)
- *optional*: `doxygen` with `dot` bindings for documentiation
The following R packages (and their dependencies) must also be
installed:
The following R libraries must then be installed, which will get the
needed dependencies automatically:
- [Rcpp](https://cran.r-project.org/web/packages/Rcpp/index.html)
- [RInside](https://cran.r-project.org/web/packages/RInside/index.html)
- [qs](https://cran.r-project.org/web/packages/qs/index.html)
- [qs2](https://cran.r-project.org/web/packages/qs2/index.html)
This can be simply achieved by issuing the following commands:
```sh
# start R environment
$ R
# install R dependencies (case sensitive!)
> install.packages(c("Rcpp", "RInside","qs","qs2"))
> q(save="no")
```
### Clone the repository
POET can be anonimously cloned from this repo over https. Make sure to
also download the submodules:
```sh
git clone --recurse-submodules https://git.gfz.de/naaice/poet.git
```
The `--recurse-submodules` option is a shorthand for:
```sh
cd poet
git submodule init && git submodule update
```
### Compiling source code
POET is built with CMake. You can generate Makefiles by running the
usual:
The generation of makefiles is done with CMake. You should be able to generate
Makefiles by running:
```sh
mkdir build && cd build
cmake -DCMAKE_BUILD_TYPE=Release ..
cmake ..
```
This will create the directory `build` and processes the CMake files
and generate Makefiles from it. You're now able to run `make` to start
build process.
If everything went well you'll find the executables at
`build/src/poet`, but it is recommended to install the POET project
If everything went well you'll find the executable at
`build/app/poet`, but it is recommended to install the POET project
structure to a desired `CMAKE_INSTALL_PREFIX` with `make install`.
During the generation of Makefiles, various options can be specified
via `cmake -D <option>=<value> [...]`. Currently, there are the
following available options:
- **POET_DHT_Debug**=_boolean_ - toggles the output of detailed
statistics about DHT usage. Defaults to _OFF_.
- **POET_ENABLE_TESTING**=_boolean_ - enables small set of unit tests
(more to come). Defaults to _OFF_.
- **POET_PHT_ADDITIONAL_INFO**=_boolean_ - enabling the count of
accesses to one PHT bucket. Use with caution, as things will get
slowed down significantly. Defaults to _OFF_.
- **POET_PREPROCESS_BENCHS**=*boolean* - enables the preprocessing of
predefined models/benchmarks. Defaults to *ON*.
- **POET_DHT_Debug**=_boolean_ - toggles the output of detailed statistics about
DHT usage. Defaults to _OFF_.
- **POET_ENABLE_TESTING**=_boolean_ - enables small set of unit tests (more to
come). Defaults to _OFF_.
- **POET_PHT_ADDITIONAL_INFO**=_boolean_ - enabling the count of accesses to one
PHT bucket. Use with caution, as things will get slowed down significantly.
Defaults to _OFF_.
- **POET_PREPROCESS_BENCHS**=*boolean* - enables the preprocessing of predefined
models/benchmarks. Defaults to *ON*.
### Example: Build from scratch
@ -115,7 +87,7 @@ follows:
$ R
# install R dependencies
> install.packages(c("Rcpp", "RInside","qs","qs2"))
> install.packages(c("Rcpp", "RInside"))
> q(save="no")
# cd into POET project root
@ -123,7 +95,7 @@ $ cd <POET_dir>
# Build process
$ mkdir build && cd build
$ cmake -DCMAKE_INSTALL_PREFIX=/home/<user>/poet -DCMAKE_BUILD_TYPE=Release ..
$ cmake -DCMAKE_INSTALL_PREFIX=/home/<user>/poet ..
$ make -j<max_numprocs>
$ make install
```
@ -143,59 +115,53 @@ poet
└── share
└── poet
├── barite
│   ├── barite_200.qs2
│   ├── barite_200.rds
│   ├── barite_200_rt.R
│   ├── barite_het.qs2
│   ├── barite_het.rds
│   └── barite_het_rt.R
├── dolo
│   ├── dolo_inner_large.qs2
│   ├── dolo_inner_large.rds
│   ├── dolo_inner_large_rt.R
│   ├── dolo_interp.qs2
│   ├── dolo_interp.rds
│   └── dolo_interp_rt.R
└── surfex
├── PoetEGU_surfex_500.qs2
├── PoetEGU_surfex_500.rds
└── PoetEGU_surfex_500_rt.R
```
With the installation of POET, two executables are provided:
- `poet` - the main executable to run simulations
- `poet_init` - a preprocessor to generate input files for POET from
R scripts
- `poet_init` - a preprocessor to generate input files for POET from R scripts
Preprocessed benchmarks can be found in the `share/poet` directory
with an according *runtime* setup. More on those files and how to
create them later.
Preprocessed benchmarks can be found in the `share/poet` directory with an
according *runtime* setup. More on those files and how to create them later.
## Running
Run POET by `mpirun ./poet [OPTIONS] <RUNFILE> <SIMFILE>
<OUTPUT_DIRECTORY>` where:
Run POET by `mpirun ./poet [OPTIONS] <RUNFILE> <SIMFILE> <OUTPUT_DIRECTORY>`
where:
- **OPTIONS** - POET options (explained below)
- **RUNFILE** - Runtime parameters described as R script
- **RUNFILE** - Runtime parameters described as R script
- **SIMFILE** - Simulation input prepared by `poet_init`
- **OUTPUT_DIRECTORY** - path, where all output of POET should be
stored
- **OUTPUT_DIRECTORY** - path, where all output of POET should be stored
### POET command line arguments
### POET options
The following parameters can be set:
| Option | Value | Description |
|-----------------------------|--------------|----------------------------------------------------------------------------------|
| **--work-package-size=** | _1..n_ | size of work packages (defaults to _5_) |
| **-P, --progress** | | show progress bar |
| **--ai-surrogate** | | activates the AI surrogate chemistry model (defaults to _OFF_) |
| **--dht** | | enabling DHT usage (defaults to _OFF_) |
| **--qs** | | store results using qs::qsave() (.qs extension) instead of default qs2 (.qs2) |
| **--rds** | | store results using saveRDS() (.rds extension) instead of default qs2 (.qs2) |
| **--dht-strategy=** | _0-1_ | change DHT strategy. **NOT IMPLEMENTED YET** (Defaults to _0_) |
| **--dht-size=** | _1-n_ | size of DHT per process involved in megabyte (defaults to _1000 MByte_) |
| **--dht-snaps=** | _0-2_ | disable or enable storage of DHT snapshots |
| **--dht-file=** | `<SNAPSHOT>` | initializes DHT with the given snapshot file |
| **--interp-size** | _1-n_ | size of PHT (interpolation) per process in megabyte |
| **--interp-bucket-entries** | _1-n_ | number of entries to store at maximum in one PHT bucket |
| **--interp-min** | _1-n_ | number of entries in PHT bucket needed to start interpolation |
| Option | Value | Description |
|-----------------------------|--------------|--------------------------------------------------------------------------------------------------------------------------|
| **--work-package-size=** | _1..n_ | size of work packages (defaults to _5_) |
| **-P, --progress** | | show progress bar |
| **--dht** | | enabling DHT usage (defaults to _OFF_) |
| **--dht-strategy=** | _0-1_ | change DHT strategy. **NOT IMPLEMENTED YET** (Defaults to _0_) |
| **--dht-size=** | _1-n_ | size of DHT per process involved in megabyte (defaults to _1000 MByte_) |
| **--dht-snaps=** | _0-2_ | disable or enable storage of DHT snapshots |
| **--dht-file=** | `<SNAPSHOT>` | initializes DHT with the given snapshot file |
| **--interp-size** | _1-n_ | size of PHT (interpolation) per process in megabyte |
| **--interp-bucket-entries** | _1-n_ | number of entries to store at maximum in one PHT bucket |
| **--interp-min** | _1-n_ | number of entries in PHT bucket needed to start interpolation |
#### Additions to `dht-snaps`
@ -209,91 +175,39 @@ Following values can be set:
### Example: Running from scratch
We will continue the above example and start a simulation with
*barite_het*, which simulation files can be found in
`<INSTALL_DIR>/share/poet/barite/barite_het*`. As transport a
heterogeneous diffusion is used. It's a small 2D grid, 2x5 grid,
simulating 50 time steps with a time step size of 100 seconds. To
start the simulation with 4 processes `cd` into your previously
installed POET-dir `<POET_INSTALL_DIR>/bin` and run:
We will continue the above example and start a simulation with *barite_het*,
which simulation files can be found in
`<INSTALL_DIR>/share/poet/barite/barite_het*`. As transport a heterogeneous
diffusion is used. It's a small 2D grid, 2x5 grid, simulating 50 time steps with
a time step size of 100 seconds. To start the simulation with 4 processes `cd`
into your previously installed POET-dir `<POET_INSTALL_DIR>/bin` and run:
```sh
cp ../share/poet/barite/barite_het* .
mpirun -n 4 ./poet barite_het_rt.R barite_het.qs2 output
mpirun -n 4 ./poet barite_het_rt.R barite_het.rds output
```
After a finished simulation all data generated by POET will be found
in the directory `output`.
You might want to use the DHT to cache previously simulated data and
reuse them in further time-steps. Just append `--dht` to the options
of POET to activate the usage of the DHT. Also, after each iteration a
DHT snapshot shall be produced. This is done by appending the
`--dht-snaps=<value>` option. The resulting call would look like this:
You might want to use the DHT to cache previously simulated data and reuse them
in further time-steps. Just append `--dht` to the options of POET to activate
the usage of the DHT. Also, after each iteration a DHT snapshot shall be
produced. This is done by appending the `--dht-snaps=<value>` option. The
resulting call would look like this:
```sh
mpirun -n 4 ./poet --dht --dht-snaps=2 barite_het_rt.R barite_het.qs2 output
mpirun -n 4 ./poet --dht --dht-snaps=2 barite_het_rt.R barite_het.rds output
```
### Example: Preparing Environment and Running with AI surrogate
To run the AI surrogate, you need to install the R package `keras3`. The
compilation process of POET remains the same as shown above.
In the following code block, the installation process on the Turing Cluster is
shown. `miniconda` is used to create a virtual environment to install
tensorflow/keras. Please adapt the installation process to your needs.
<!-- Start an R interactive session and install the required packages: -->
```sh
# First, install the required R packages
R -e "install.packages('keras3', repos='https://cloud.r-project.org/')"
# manually create a virtual environment to install keras/python using conda,
# as this is somehow broken on the Turing Cluster when using the `keras::install_keras()` function
cd poet
# create a virtual environment in the .ai directory with python 3.11
conda create -p ./.ai python=3.11
conda activate ./.ai
# install tensorflow and keras
pip install keras tensorflow[and-cuda]
# add conda's python path to the R environment
# make sure to have the conda environment activated
echo -e "RETICULATE_PYTHON=$(which python)\n" >> ~/.Renviron
```
After setup the R environment, recompile POET and you're ready to run the AI
surrogate.
```sh
cd <installation_dir>/bin
# copy the benchmark files to the installation directory
cp <project_root_dir>/bench/barite/{barite_50ai*,db_barite.dat,barite.pqi} .
# preprocess the benchmark
./poet_init barite_50ai.R
# run POET with AI surrogate and GPU utilization
srun --gres=gpu -N 1 -n 12 ./poet --ai-surrogate barite_50ai_rt.R barite_50ai.qs2 output
```
Keep in mind that the AI surrogate is currently not stable or might also not
produce any valid predictions.
## Defining a model
In order to provide a model to POET, you need to setup a R script
which can then be used by `poet_init` to generate the simulation
input. Which parameters are required can be found in the
[Wiki](https://git.gfz.de/naaice/poet/-/wikis/Initialization).
We try to keep the document up-to-date. However, if you encounter
missing information or need help, please get in touch with us via the
issue tracker or E-Mail.
In order to provide a model to POET, you need to setup a R script which can then
be used by `poet_init` to generate the simulation input. Which parameters are
required can be found in the
[Wiki](https://git.gfz-potsdam.de/naaice/poet/-/wikis/Initialization). We try to
keep the document up-to-date. However, if you encounter missing information or
need help, please get in touch with us via the issue tracker or E-Mail.
`poet_init` can be used as follows:
@ -303,12 +217,11 @@ issue tracker or E-Mail.
where:
- **output** - name of the output file (defaults to the input file
name with the extension `.qs2`)
- **setwd** - set the working directory to the directory of the input
file (e.g. to allow relative paths in the input script). However,
the output file will be stored in the directory from which
`poet_init` was called.
- **output** - name of the output file (defaults to the input file name
with the extension `.rds`)
- **setwd** - set the working directory to the directory of the input file (e.g.
to allow relative paths in the input script). However, the output file
will be stored in the directory from which `poet_init` was called.
## About the usage of MPI_Wtime()
@ -318,44 +231,3 @@ important information from the OpenMPI Man Page:
For example, on platforms that support it, the clock_gettime()
function will be used to obtain a monotonic clock value with whatever
precision is supported on that platform (e.g., nanoseconds).
## Additional functions for the AI surrogate
The AI surrogate can be activated for any benchmark and is by default
initiated as a sequential keras model with three hidden layer of depth
48, 96, 24 with relu activation and adam optimizer. All functions in
`ai_surrogate_model.R` can be overridden by adding custom definitions
via an R file in the input script. This is done by adding the path to
this file in the input script. Simply add the path as an element
called `ai_surrogate_input_script` to the `chemistry_setup` list.
Please use the global variable `ai_surrogate_base_path` as a base path
when relative filepaths are used in custom funtions.
**There is currently no default implementation to determine the
validity of predicted values.** This means, that every input script
must include an R source file with a custom function
`validate_predictions(predictors, prediction)`. Examples for custom
functions can be found for the barite_200 benchmark
The functions can be defined as follows:
`validate_predictions(predictors, prediction)`: Returns a boolean
index vector that signals for each row in the predictions if the
values are considered valid. Can eg. be implemented as a mass balance
threshold between the predictors and the prediction.
`initiate_model()`: Returns a keras model. Can be used to load
pretrained models.
`preprocess(df, backtransform = FALSE, outputs = FALSE)`: Returns the
scaled/transformed/backtransformed dataframe. The `backtransform` flag
signals if the current processing step is applied to data that's
assumed to be scaled and expects backtransformed values. The `outputs`
flag signals if the current processing step is applied to the output
or tatget of the model. This can be used to eg. skip these processing
steps and only scale the model input.
`training_step (model, predictor, target, validity)`: Trains the model
after each iteration. `validity` is the bool index vector given by
`validate_predictions` and can eg. be used to only train on values
that have not been valid predictions.

75
R_lib/ai_surrogate_model.R
View File

@ -1,75 +0,0 @@
## This file contains default function implementations for the ai surrogate.
## To load pretrained models, use pre-/postprocessing or change hyperparameters
## it is recommended to override these functions with custom implementations via
## the input script. The path to the R-file containing the functions mus be set
## in the variable "ai_surrogate_input_script". See the barite_200.R file as an
## example and the general README for more information.
require(keras3)
require(tensorflow)
initiate_model <- function() {
hidden_layers <- c(48, 96, 24)
activation <- "relu"
loss <- "mean_squared_error"
input_length <- length(ai_surrogate_species)
output_length <- length(ai_surrogate_species)
## Creates a new sequential model from scratch
model <- keras_model_sequential()
## Input layer defined by input data shape
model %>% layer_dense(units = input_length,
activation = activation,
input_shape = input_length,
dtype = "float32")
for (layer_size in hidden_layers) {
model %>% layer_dense(units = layer_size,
activation = activation,
dtype = "float32")
}
## Output data defined by output data shape
model %>% layer_dense(units = output_length,
activation = activation,
dtype = "float32")
model %>% compile(loss = loss,
optimizer = "adam")
return(model)
}
gpu_info <- function() {
msgm(tf_gpu_configured())
}
prediction_step <- function(model, predictors) {
prediction <- predict(model, as.matrix(predictors))
colnames(prediction) <- colnames(predictors)
return(as.data.frame(prediction))
}
preprocess <- function(df, backtransform = FALSE, outputs = FALSE) {
return(df)
}
postprocess <- function(df, backtransform = TRUE, outputs = TRUE) {
return(df)
}
set_valid_predictions <- function(temp_field, prediction, validity) {
temp_field[validity == 1, ] <- prediction[validity == 1, ]
return(temp_field)
}
training_step <- function(model, predictor, target, validity) {
msgm("Training:")
x <- as.matrix(predictor)
y <- as.matrix(target[colnames(x)])
model %>% fit(x, y)
model %>% save_model_tf(paste0(out_dir, "/current_model.keras"))
}

72
R_lib/init_r_lib.R
View File

@ -1,67 +1,21 @@
### Copyright (C) 2018-2024 Marco De Lucia, Max Luebke (GFZ Potsdam, University of Potsdam)
###
### POET is free software; you can redistribute it and/or modify it under the
### terms of the GNU General Public License as published by the Free Software
### Foundation; either version 2 of the License, or (at your option) any later
### version.
###
### POET is distributed in the hope that it will be useful, but WITHOUT ANY
### WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
### A PARTICULAR PURPOSE. See the GNU General Public License for more details.
###
### You should have received a copy of the GNU General Public License along with
### this program; if not, write to the Free Software Foundation, Inc., 51
### Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
##' @param pqc_mat matrix, containing IDs and PHREEQC outputs
##' @param grid matrix, zonation referring to pqc_mat$ID
##' @return a data.frame
# pqc_to_grid <- function(pqc_mat, grid) {
# # Convert the input DataFrame to a matrix
# pqc_mat <- as.matrix(pqc_mat)
# # Flatten the matrix into a vector
# id_vector <- as.integer(t(grid))
# # Find the matching rows in the matrix
# row_indices <- match(id_vector, pqc_mat[, "ID"])
# # Extract the matching rows from pqc_mat to size of grid matrix
# result_mat <- pqc_mat[row_indices, ]
# # Convert the result matrix to a data frame
# res_df <- as.data.frame(result_mat)
# # Remove all columns which only contain NaN
# res_df <- res_df[, colSums(is.na(res_df)) != nrow(res_df)]
# # Remove row names
# rownames(res_df) <- NULL
# return(res_df)
# }
##' @param pqc_mat matrix, containing IDs and PHREEQC outputs
##' @param grid matrix, zonation referring to pqc_mat$ID
##' @return a data.frame
pqc_to_grid <- function(pqc_mat, grid) {
pqc_to_grid <- function(pqc_in, grid) {
# Convert the input DataFrame to a matrix
pqc_mat <- as.matrix(pqc_mat)
pqc_in <- as.matrix(pqc_in)
# Flatten the matrix into a vector
id_vector <- as.integer(t(grid))
id_vector <- as.numeric(t(grid))
# Find the matching rows in the matrix
row_indices <- match(id_vector, pqc_mat[, "ID"])
row_indices <- match(id_vector, pqc_in[, "ID"])
# Extract the matching rows from pqc_mat to size of grid matrix
result_mat <- pqc_mat[row_indices, ]
# Extract the matching rows from pqc_in to size of grid matrix
result_mat <- pqc_in[row_indices, ]
# Convert the result matrix to a data frame
res_df <- as.data.frame(result_mat)
# Remove all columns which only contain NaN
# res_df <- res_df[, colSums(is.na(res_df)) != nrow(res_df)]
res_df <- res_df[, colSums(is.na(res_df)) != nrow(res_df)]
# Remove row names
rownames(res_df) <- NULL
@ -69,12 +23,6 @@ pqc_to_grid <- function(pqc_mat, grid) {
return(res_df)
}
##' @param pqc_mat matrix,
##' @param transport_spec column name of species in pqc_mat
##' @param id
##' @title
##' @return
resolve_pqc_bound <- function(pqc_mat, transport_spec, id) {
df <- as.data.frame(pqc_mat, check.names = FALSE)
value <- df[df$ID == id, transport_spec]
@ -86,10 +34,6 @@ resolve_pqc_bound <- function(pqc_mat, transport_spec, id) {
return(value)
}
##' @title
##' @param init_grid
##' @param new_names
##' @return
add_missing_transport_species <- function(init_grid, new_names) {
# add 'ID' to new_names front, as it is not a transport species but required
new_names <- c("ID", new_names)
@ -109,4 +53,4 @@ add_missing_transport_species <- function(init_grid, new_names) {
new_grid <- cbind(new_grid, append_df)
return(new_grid)
}
}

190
R_lib/kin_r_library.R
View File

@ -1,3 +1,5 @@
## Time-stamp: "Last modified 2023-08-15 11:58:23 delucia"
### Copyright (C) 2018-2023 Marco De Lucia, Max Luebke (GFZ Potsdam)
###
### POET is free software; you can redistribute it and/or modify it under the
@ -33,18 +35,14 @@ master_init <- function(setup, out_dir, init_field) {
setup$iterations <- setup$maxiter
setup$simulation_time <- 0
dgts <- as.integer(ceiling(log10(setup$maxiter)))
## string format to use in sprintf
fmt <- paste0("%0", dgts, "d")
if (is.null(setup[["store_result"]])) {
setup$store_result <- TRUE
}
if (setup$store_result) {
init_field_out <- paste0(out_dir, "/iter_", sprintf(fmt = fmt, 0), ".", setup$out_ext)
init_field_out <- paste0(out_dir, "/iter_0.rds")
init_field <- data.frame(init_field, check.names = FALSE)
SaveRObj(x = init_field, path = init_field_out)
saveRDS(init_field, file = init_field_out)
msgm("Stored initial field in ", init_field_out)
if (is.null(setup[["out_save"]])) {
setup$out_save <- seq(1, setup$iterations)
@ -63,7 +61,7 @@ master_iteration_end <- function(setup, state_T, state_C) {
iter <- setup$iter
# print(iter)
## max digits for iterations
dgts <- as.integer(ceiling(log10(setup$maxiter + 1)))
dgts <- as.integer(ceiling(log10(setup$maxiter)))
## string format to use in sprintf
fmt <- paste0("%0", dgts, "d")
@ -71,30 +69,23 @@ master_iteration_end <- function(setup, state_T, state_C) {
## comprised in setup$out_save
if (setup$store_result) {
if (iter %in% setup$out_save) {
nameout <- paste0(setup$out_dir, "/iter_", sprintf(fmt = fmt, iter), ".", setup$out_ext)
nameout <- paste0(setup$out_dir, "/iter_", sprintf(fmt = fmt, iter), ".rds")
state_T <- data.frame(state_T, check.names = FALSE)
state_C <- data.frame(state_C, check.names = FALSE)
ai_surrogate_info <- list(
prediction_time = if (exists("ai_prediction_time")) as.integer(ai_prediction_time) else NULL,
training_time = if (exists("ai_training_time")) as.integer(ai_training_time) else NULL,
valid_predictions = if (exists("validity_vector")) validity_vector else NULL
)
SaveRObj(x = list(
saveRDS(list(
T = state_T,
C = state_C,
simtime = as.integer(setup$simulation_time),
totaltime = as.integer(totaltime),
ai_surrogate_info = ai_surrogate_info
), path = nameout)
simtime = as.integer(setup$simulation_time)
), file = nameout)
msgm("results stored in <", nameout, ">")
}
}
## Add last time step to simulation time
setup$simulation_time <- setup$simulation_time + setup$timesteps[iter]
## msgm("done iteration", iter, "/", length(setup$timesteps))
msgm("done iteration", iter, "/", length(setup$timesteps))
setup$iter <- setup$iter + 1
return(setup)
}
@ -109,103 +100,70 @@ msgm <- function(...) {
}
## Function called by master R process to store on disk all relevant
## parameters for the simulation
StoreSetup <- function(setup, filesim, out_dir) {
to_store <- vector(mode = "list", length = 4)
## names(to_store) <- c("Sim", "Flow", "Transport", "Chemistry", "DHT")
names(to_store) <- c("Sim", "Transport", "DHT", "Cmdline")
## read the setup R file, which is sourced in kin.cpp
tmpbuff <- file(filesim, "r")
setupfile <- readLines(tmpbuff)
close.connection(tmpbuff)
to_store$Sim <- setupfile
## to_store$Flow <- list(
## snapshots = setup$snapshots,
## gridfile = setup$gridfile,
## phase = setup$phase,
## density = setup$density,
## dt_differ = setup$dt_differ,
## prolong = setup$prolong,
## maxiter = setup$maxiter,
## saved_iter = setup$iter_output,
## out_save = setup$out_save )
to_store$Transport <- setup$diffusion
## to_store$Chemistry <- list(
## nprocs = n_procs,
## wp_size = work_package_size,
## base = setup$base,
## first = setup$first,
## init = setup$initsim,
## db = db,
## kin = setup$kin,
## ann = setup$ann)
if (dht_enabled) {
to_store$DHT <- list(
enabled = dht_enabled,
log = dht_log
## signif = dht_final_signif,
## proptype = dht_final_proptype
)
} else {
to_store$DHT <- FALSE
}
if (dht_enabled) {
to_store$DHT <- list(
enabled = dht_enabled,
log = dht_log
# signif = dht_final_signif,
# proptype = dht_final_proptype
)
} else {
to_store$DHT <- FALSE
}
saveRDS(to_store, file = paste0(fileout, "/setup.rds"))
msgm("initialization stored in ", paste0(fileout, "/setup.rds"))
}
GetWorkPackageSizesVector <- function(n_packages, package_size, len) {
ids <- rep(1:n_packages, times = package_size, each = 1)[1:len]
return(as.integer(table(ids)))
}
## Handler to read R objs from binary files using either builtin
## readRDS(), qs::qread() or qs2::qs_read() based on file extension
ReadRObj <- function(path) {
## code borrowed from tools::file_ext()
pos <- regexpr("\\.([[:alnum:]]+)$", path)
extension <- ifelse(pos > -1L, substring(path, pos + 1L), "")
switch(extension,
rds = readRDS(path),
qs = qs::qread(path),
qs2 = qs2::qs_read(path)
)
}
## Handler to store R objs to binary files using either builtin
## saveRDS() or qs::qsave() based on file extension
SaveRObj <- function(x, path) {
## msgm("Storing to", path)
## code borrowed from tools::file_ext()
pos <- regexpr("\\.([[:alnum:]]+)$", path)
extension <- ifelse(pos > -1L, substring(path, pos + 1L), "")
switch(extension,
rds = saveRDS(object = x, file = path),
qs = qs::qsave(x = x, file = path),
qs2 = qs2::qs_save(object = x, file = path)
)
}
######## Old relic code
## ## Function called by master R process to store on disk all relevant
## ## parameters for the simulation
## StoreSetup <- function(setup, filesim, out_dir) {
## to_store <- vector(mode = "list", length = 4)
## ## names(to_store) <- c("Sim", "Flow", "Transport", "Chemistry", "DHT")
## names(to_store) <- c("Sim", "Transport", "DHT", "Cmdline")
## ## read the setup R file, which is sourced in kin.cpp
## tmpbuff <- file(filesim, "r")
## setupfile <- readLines(tmpbuff)
## close.connection(tmpbuff)
## to_store$Sim <- setupfile
## ## to_store$Flow <- list(
## ## snapshots = setup$snapshots,
## ## gridfile = setup$gridfile,
## ## phase = setup$phase,
## ## density = setup$density,
## ## dt_differ = setup$dt_differ,
## ## prolong = setup$prolong,
## ## maxiter = setup$maxiter,
## ## saved_iter = setup$iter_output,
## ## out_save = setup$out_save )
## to_store$Transport <- setup$diffusion
## ## to_store$Chemistry <- list(
## ## nprocs = n_procs,
## ## wp_size = work_package_size,
## ## base = setup$base,
## ## first = setup$first,
## ## init = setup$initsim,
## ## db = db,
## ## kin = setup$kin,
## ## ann = setup$ann)
## if (dht_enabled) {
## to_store$DHT <- list(
## enabled = dht_enabled,
## log = dht_log
## ## signif = dht_final_signif,
## ## proptype = dht_final_proptype
## )
## } else {
## to_store$DHT <- FALSE
## }
## if (dht_enabled) {
## to_store$DHT <- list(
## enabled = dht_enabled,
## log = dht_log
## # signif = dht_final_signif,
## # proptype = dht_final_proptype
## )
## } else {
## to_store$DHT <- FALSE
## }
## saveRDS(to_store, file = paste0(fileout, "/setup.rds"))
## msgm("initialization stored in ", paste0(fileout, "/setup.rds"))
## }

11
bench/CMakeLists.txt
View File

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

9
bench/barite/README.org
View File

@ -22,15 +22,6 @@ mpirun -np 4 ./poet --interp barite_interp_eval.R barite_results
grid
- =barite_200.R=: POET input script for a 200x200 simulation
grid
- =barite_200ai_surrogate_input_script.R=: Defines the ai surrogate functions
to load a pretrained model and apply min-max-feature scaling on the model inputs
and target. Prediction validity is assessed with a threshold of 3e-5 on the mass
balance of Ba and Sr.
- =barite_200min_max_bounds=: Minimum and maximum values from 50 iterations of the
barite_200 benchmark. Used for feature scaling in the ai surrogate.
- =barite_200model_min_max.keras=: A sequential keras model that has been trained
on 50 iterations of the barite_200 benchmark with min-max-scaled inputs
and targets/outputs.
- =db_barite.dat=: PHREEQC database containing the kinetic expressions
for barite and celestite, stripped down from =phreeqc.dat=
- =barite.pqi=: PHREEQC input script defining the chemical system

21
bench/barite/barite_200.R
View File

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

48
bench/barite/barite_200ai_surrogate_input_script.R
View File

@ -1,48 +0,0 @@
## load a pretrained model from tensorflow file
## Use the global variable "ai_surrogate_base_path" when using file paths
## relative to the input script
initiate_model <- function() {
init_model <- normalizePath(paste0(ai_surrogate_base_path,
"model_min_max_float64.keras"))
return(load_model_tf(init_model))
}
scale_min_max <- function(x, min, max, backtransform) {
if (backtransform) {
return((x * (max - min)) + min)
} else {
return((x - min) / (max - min))
}
}
preprocess <- function(df, backtransform = FALSE, outputs = FALSE) {
minmax_file <- normalizePath(paste0(ai_surrogate_base_path,
"min_max_bounds.rds"))
global_minmax <- readRDS(minmax_file)
for (column in colnames(df)) {
df[column] <- lapply(df[column],
scale_min_max,
global_minmax$min[column],
global_minmax$max[column],
backtransform)
}
return(df)
}
mass_balance <- function(predictors, prediction) {
dBa <- abs(prediction$Ba + prediction$Barite -
predictors$Ba - predictors$Barite)
dSr <- abs(prediction$Sr + prediction$Celestite -
predictors$Sr - predictors$Celestite)
return(dBa + dSr)
}
validate_predictions <- function(predictors, prediction) {
epsilon <- 3e-5
mb <- mass_balance(predictors, prediction)
msgm("Mass balance mean:", mean(mb))
msgm("Mass balance variance:", var(mb))
msgm("Rows where mass balance meets threshold", epsilon, ":",
sum(mb < epsilon))
return(mb < epsilon)
}

60
bench/barite/barite_50ai.R
View File

@ -1,60 +0,0 @@
## Time-stamp: "Last modified 2024-05-30 13:34:14 delucia"
cols <- 50
rows <- 50
s_cols <- 0.25
s_rows <- 0.25
grid_def <- matrix(2, nrow = rows, ncol = cols)
# Define grid configuration for POET model
grid_setup <- list(
pqc_in_file = "./barite.pqi",
pqc_db_file = "./db_barite.dat", ## Path to the database file for Phreeqc
grid_def = grid_def, ## Definition of the grid, containing IDs according to the Phreeqc input script
grid_size = c(s_rows, s_cols), ## Size of the grid in meters
constant_cells = c() ## IDs of cells with constant concentration
)
bound_length <- 2
bound_def <- list(
"type" = rep("constant", bound_length),
"sol_id" = rep(3, bound_length),
"cell" = seq(1, bound_length)
)
homogenous_alpha <- 1e-8
diffusion_setup <- list(
boundaries = list(
"W" = bound_def,
"N" = bound_def
),
alpha_x = homogenous_alpha,
alpha_y = homogenous_alpha
)
dht_species <- c(
"H" = 3,
"O" = 3,
"Charge" = 3,
"Ba" = 6,
"Cl" = 6,
"S" = 6,
"Sr" = 6,
"Barite" = 5,
"Celestite" = 5
)
chemistry_setup <- list(
dht_species = dht_species,
ai_surrogate_input_script = "./barite_50ai_surr_mdl.R"
)
# Define a setup list for simulation configuration
setup <- list(
Grid = grid_setup, # Parameters related to the grid structure
Diffusion = diffusion_setup, # Parameters related to the diffusion process
Chemistry = chemistry_setup
)

BIN
bench/barite/barite_50ai_all.keras
View File

Binary file not shown.

9
bench/barite/barite_50ai_rt.R
View File

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

90
bench/barite/barite_50ai_surr_mdl.R
View File

@ -1,90 +0,0 @@
## Time-stamp: "Last modified 2024-05-30 13:27:06 delucia"
## load a pretrained model from tensorflow file
## Use the global variable "ai_surrogate_base_path" when using file paths
## relative to the input script
initiate_model <- function() {
require(keras3)
require(tensorflow)
init_model <- normalizePath(paste0(ai_surrogate_base_path,
"barite_50ai_all.keras"))
Model <- keras3::load_model(init_model)
msgm("Loaded model:")
print(str(Model))
return(Model)
}
scale_min_max <- function(x, min, max, backtransform) {
if (backtransform) {
return((x * (max - min)) + min)
} else {
return((x - min) / (max - min))
}
}
minmax <- list(min = c(H = 111.012433592824, O = 55.5062185549492, Charge = -3.1028354471876e-08,
Ba = 1.87312878574393e-141, Cl = 0, `S(6)` = 4.24227510643685e-07,
Sr = 0.00049382996130541, Barite = 0.000999542409828586, Celestite = 0.244801877115968),
max = c(H = 111.012433679682, O = 55.5087003521685, Charge = 5.27666636082035e-07,
Ba = 0.0908849779513762, Cl = 0.195697626449355, `S(6)` = 0.000620774752665846,
Sr = 0.0558680070692722, Barite = 0.756779139057097, Celestite = 1.00075422160624
))
preprocess <- function(df) {
if (!is.data.frame(df))
df <- as.data.frame(df, check.names = FALSE)
as.data.frame(lapply(colnames(df),
function(x) scale_min_max(x=df[x],
min=minmax$min[x],
max=minmax$max[x],
backtransform=FALSE)),
check.names = FALSE)
}
postprocess <- function(df) {
if (!is.data.frame(df))
df <- as.data.frame(df, check.names = FALSE)
as.data.frame(lapply(colnames(df),
function(x) scale_min_max(x=df[x],
min=minmax$min[x],
max=minmax$max[x],
backtransform=TRUE)),
check.names = FALSE)
}
mass_balance <- function(predictors, prediction) {
dBa <- abs(prediction$Ba + prediction$Barite -
predictors$Ba - predictors$Barite)
dSr <- abs(prediction$Sr + prediction$Celestite -
predictors$Sr - predictors$Celestite)
return(dBa + dSr)
}
validate_predictions <- function(predictors, prediction) {
epsilon <- 1E-7
mb <- mass_balance(predictors, prediction)
msgm("Mass balance mean:", mean(mb))
msgm("Mass balance variance:", var(mb))
ret <- mb < epsilon
msgm("Rows where mass balance meets threshold", epsilon, ":",
sum(ret))
return(ret)
}
training_step <- function(model, predictor, target, validity) {
msgm("Starting incremental training:")
## x <- as.matrix(predictor)
## y <- as.matrix(target[colnames(x)])
history <- model %>% keras3::fit(x = data.matrix(predictor),
y = data.matrix(target),
epochs = 10, verbose=1)
keras3::save_model(model,
filepath = paste0(out_dir, "/current_model.keras"),
overwrite=TRUE)
return(model)
}

BIN
bench/barite/min_max_bounds.rds
View File

Binary file not shown.

BIN
bench/barite/model_min_max_float64.keras
View File

Binary file not shown.

48
bench/dolo/dolo_interp.R
View File

@ -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" = 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(4)" = 6,
"Ca" = 6,
"Cl" = 3,
"Mg" = 5,
"Calcite" = 4,
"Dolomite" = 4
)
return(input[names(dht_species)])
}
check_sign_cal_dol_interp <- function(to_interp, data_set) {
dht_species <- c(
"H" = 3,
"O" = 3,
"Charge" = 3,
"C" = 6,
"C(4)" = 6,
"Ca" = 6,
"Cl" = 3,
"Mg" = 5,
@ -95,7 +95,7 @@ check_neg_cal_dol <- function(result) {
# significant digits)
pht_species <- c(
"C" = 3,
"C(4)" = 3,
"Ca" = 3,
"Mg" = 2,
"Calcite" = 2,
@ -107,7 +107,7 @@ chemistry_setup <- list(
"H" = 3,
"O" = 3,
"Charge" = 3,
"C" = 6,
"C(4)" = 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
)

BIN
bench/fgcs/20241211_README.pdf
View File

Binary file not shown.

102
bench/fgcs/20241211_README.tex
View File

@ -1,102 +0,0 @@
% Created 2024-12-11 mer 23:24
% Intended LaTeX compiler: pdflatex
\documentclass[a4paper, 9pt]{article}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{graphicx}
\usepackage{longtable}
\usepackage{wrapfig}
\usepackage{rotating}
\usepackage[normalem]{ulem}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{capt-of}
\usepackage{hyperref}
\usepackage{fullpage}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{charter}
\usepackage{listings}
\lstloadlanguages{R}
\author{MDL <delucia@gfz.de>}
\date{2024-12-11}
\title{A \texttt{barite}-based benchmark for FGCS interpolation paper}
\begin{document}
\maketitle
\section{Description}
\label{sec:org739879a}
\begin{itemize}
\item \texttt{barite\_fgcs\_2.R}: POET input script with circular
"crystals" on a 200x200 nodes grid
\item \(\alpha\): isotropic 10\textsuperscript{-5}
m\textsuperscript{2}/s outside of the crystals,
10\textsuperscript{-7} inside
\item 200 iterations, dt = 1000
\item \texttt{barite\_fgcs\_2.pqi}: PHREEQC input, 4 SOLUTIONS
(basically the same as in \texttt{barite} benchmark):
\begin{enumerate}
\item Equilibrium with Celestite, no mineral \(Rightarrow\)
\item Equilibrium with Celestite, KINETICS Celestite (1 mol) and
Barite (0 mol)
\item Injection of 0.1 BaCl2 from NW corner
\item Injection of 0.2 BaCl2 from SE corner
\end{enumerate}
\item \texttt{db\_barite.dat}: PHREEQC database containing the kinetic
expressions for barite and celestite, stripped down from
\texttt{phreeqc.dat}
\end{itemize}
\begin{figure}[htbp]
\centering
\includegraphics[width=0.48\textwidth]{./fgcs_Celestite_init.pdf}
\includegraphics[width=0.48\textwidth]{./fgcs_Barite_200.pdf}
\caption{\textbf{Left:} Initial distribution of Celestite
"crystals". \textbf{Right:} precipitated Barite}
\end{figure}
\section{Interpolation}
\label{sec:org2a09431}
Using the following parametrization:
\begin{lstlisting}
dht_species <- c("H" = 7,
"O" = 7,
"Ba" = 7,
"Cl" = 7,
"S(6)" = 7,
"Sr" = 7,
"Barite" = 4,
"Celestite" = 4)
pht_species <- c("Ba" = 4,
"Cl" = 3,
"S(6)" = 3,
"Sr" = 3,
"Barite" = 2,
"Celestite" = 2 )
\end{lstlisting}
Runtime goes from 1800 to 600 s (21 CPUs) but there are "suspect"
errors especially in O and H, where "suspect" means some values appear
to be multiplied by 2:
\begin{figure}[htbp]
\centering
\includegraphics[width=0.9\textwidth]{./fgcs_interp_1.pdf}
\caption{Scatterplots reference vs interpolated after 1 coupling
iteration}
\end{figure}
\end{document}
%%% Local Variables:
%%% mode: LaTeX
%%% TeX-master: t
%%% End:

90
bench/fgcs/EvalFGCS.R
View File

@ -1,90 +0,0 @@
## Time-stamp: "Last modified 2024-12-11 23:21:25 delucia"
library(PoetUtils)
library(viridis)
res <- ReadPOETSims("./res_fgcs2_96/")
pp <- PlotField(res$iter_200$C$Barite, rows = 200, cols = 200, contour = FALSE,
nlevels=12, palette=terrain.colors)
cairo_pdf("fgcs_Celestite_init.pdf", family="serif")
par(mar=c(0,0,0,0))
pp <- PlotField((res$iter_000$Celestite), rows = 200, cols = 200,
contour = FALSE, breaks=c(-0.5,0.5,1.5),
palette = grey.colors, plot.axes = FALSE, scale = FALSE,
main="Initial Celestite crystals")
dev.off()
cairo_pdf("fgcs_Ba_init.pdf", family="serif")
par(mar=c(0,0,0,0))
pp <- PlotField(log10(res$iter_001$C$Cl), rows = 200, cols = 200,
contour = FALSE,
palette = terrain.colors, plot.axes = FALSE, scale = FALSE,
main="log10(Ba)")
dev.off()
pp <- PlotField(log10(res$iter_002$C$Ba), rows = 200, cols = 200,
contour = FALSE, palette = viridis, rev.palette = FALSE,
main = "log10(Ba) after 5 iterations")
pp <- PlotField(log10(res$iter_200$C$`S(6)`), rows = 200, cols = 200, contour = FALSE)
str(res$iter_00)
res$iter_178$C$Barite
pp <- res$iter_043$C$Barite
breaks <- pretty(pp, n = 5)
br <- c(0, 0.0005, 0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1)
pp <- PlotField(res$iter_200$C$Barite, rows = 200, cols = 200, contour = FALSE,
breaks = br, palette=terrain.colors)
cairo_pdf("fgcs_Barite_200.pdf", family="serif")
pp <- PlotField(log10(res$iter_200$C$Barite), rows = 200, cols = 200,
contour = FALSE, palette = terrain.colors, plot.axes = FALSE,
rev.palette = FALSE, main = "log10(Barite) after 200 iter")
dev.off()
ref <- ReadPOETSims("./res_fgcs_2_ref")
rei <- ReadPOETSims("./res_fgcs_2_interp1/")
timref <- ReadRObj("./res_fgcs_2_ref/timings.qs")
timint <- ReadRObj("./res_fgcs_2_interp1/timings.qs")
timref
timint
wch <- c("H","O", "Ba", "Sr","Cl", "S(6)")
rf <- data.matrix(ref$iter_001$C[, wch])
r1 <- data.matrix(rei$iter_001$C[, wch])
r1[is.nan(r1)] <- NA
rf[is.nan(rf)] <- NA
cairo_pdf("fgcs_interp_1.pdf", family="serif", width = 10, height = 7)
PlotScatter(rf, r1, which = wch, labs = c("ref", "interp"), cols = 3, log="", las = 1, pch=4)
dev.off()
head(r1)
head(rf)
rf$O
r1$O

2
bench/fgcs/README.org
View File

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

105
bench/fgcs/barite_fgcs_2.R
View File

@ -1,105 +0,0 @@
## Time-stamp: "Last modified 2024-12-11 16:08:11 delucia"
cols <- 1000
rows <- 1000
dim_cols <- 50
dim_rows <- 50
ncirc <- 20 ## number of crystals
rmax <- cols / 10 ## max radius (in nodes)
set.seed(22933)
centers <- cbind(sample(seq_len(cols), ncirc), sample(seq_len(rows), ncirc))
radii <- sample(seq_len(rmax), ncirc, replace = TRUE)
mi <- matrix(rep(seq_len(cols), rows), byrow = TRUE, nrow = rows)
mj <- matrix(rep(seq_len(cols), each = rows), byrow = TRUE, nrow = rows)
tmpl <- lapply(seq_len(ncirc), function(x) which((mi - centers[x, 1])^2 + (mj - centers[x, 2])^2 < radii[x]^2, arr.ind = TRUE))
inds <- do.call(rbind, tmpl)
grid <- matrix(1, nrow = rows, ncol = cols)
grid[inds] <- 2
alpha <- matrix(1e-5, ncol = cols, nrow = rows)
alpha[inds] <- 1e-7
## image(grid, asp=1)
## Define grid configuration for POET model
grid_setup <- list(
pqc_in_file = "./barite_fgcs_2.pqi",
pqc_db_file = "../barite/db_barite.dat", ## database file
grid_def = grid, ## grid definition, IDs according to the Phreeqc input
grid_size = c(dim_cols, dim_rows), ## grid size in meters
constant_cells = c() ## IDs of cells with constant concentration
)
bound_length <- cols / 10
bound_N <- list(
"type" = rep("constant", bound_length),
"sol_id" = rep(3, bound_length),
"cell" = seq(1, bound_length)
)
bound_W <- list(
"type" = rep("constant", bound_length),
"sol_id" = rep(3, bound_length),
"cell" = seq(1, bound_length)
)
bound_E <- list(
"type" = rep("constant", bound_length),
"sol_id" = rep(4, bound_length),
"cell" = seq(rows - bound_length + 1, rows)
)
bound_S <- list(
"type" = rep("constant", bound_length),
"sol_id" = rep(4, bound_length),
"cell" = seq(cols - bound_length + 1, cols)
)
diffusion_setup <- list(
boundaries = list(
"W" = bound_W,
"N" = bound_N,
"E" = bound_E,
"S" = bound_S
),
alpha_x = alpha,
alpha_y = alpha
)
dht_species <- c(
"H" = 7,
"O" = 7,
"Ba" = 7,
"Cl" = 7,
"S" = 7,
"Sr" = 7,
"Barite" = 4,
"Celestite" = 4
)
pht_species <- c(
"Ba" = 4,
"Cl" = 3,
"S" = 3,
"Sr" = 3,
"Barite" = 0,
"Celestite" = 0
)
chemistry_setup <- list(
dht_species = dht_species,
pht_species = pht_species
)
## Define a setup list for simulation configuration
setup <- list(
Grid = grid_setup, ## Parameters related to the grid structure
Diffusion = diffusion_setup, ## Parameters related to the diffusion process
Chemistry = chemistry_setup
)

49
bench/fgcs/barite_fgcs_2.pqi
View File

@ -1,49 +0,0 @@
SOLUTION 1
units mol/kgw
water 1
temperature 25
pH 7.008
pe 10.798
S 6.205e-04
Sr 6.205e-04
END
SOLUTION 2
units mol/kgw
water 1
temperature 25
pH 7.008
pe 10.798
S 6.205e-04
Sr 6.205e-04
KINETICS 2
Barite
-m 0.00
-parms 50. # reactive surface area
-tol 1e-9
Celestite
-m 1
-parms 10.0 # reactive surface area
-tol 1e-9
END
SOLUTION 3
units mol/kgw
water 1
temperature 25
Ba 0.1
Cl 0.2
END
SOLUTION 4
units mol/kgw
water 1
temperature 25
Ba 0.2
Cl 0.4
END
RUN_CELLS
-cells 1 2 3 4
END

7
bench/fgcs/barite_fgcs_2_rt.R
View File

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

637
docs/20230720_Scheme_POET_en.svg Normal file
View File

File diff suppressed because one or more lines are too long
Side by Side

After

Width:  |  Height:  |  Size: 131 KiB

281
docs/POET.drawio
View File

File diff suppressed because one or more lines are too long

4
docs/POET_scheme.svg
View File

File diff suppressed because one or more lines are too long
Side by Side

Before

Width:  |  Height:  |  Size: 813 KiB

602
docs/Scheme_POET_en.svg
View File

File diff suppressed because one or more lines are too long
Side by Side

Before

Width:  |  Height:  |  Size: 127 KiB

1
ext/iphreeqc Submodule

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

1
ext/litephreeqc

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

2
ext/tug

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

1
src/Base/Macros.hpp
View File

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

60
src/Base/RInsidePOET.hpp
View File

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

459
src/Base/argh.hpp Normal file
View File

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

30
src/CMakeLists.txt
View File

@ -1,4 +1,7 @@
add_library(POETLib
add_library(POETLib)
target_sources(POETLib
PRIVATE
Init/InitialList.cpp
Init/GridInit.cpp
Init/DiffusionInit.cpp
@ -15,15 +18,6 @@ add_library(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
@ -33,16 +27,6 @@ target_link_libraries(
PUBLIC MPI::MPI_C
)
include(FetchContent)
FetchContent_Declare(
cli11
QUIET
GIT_REPOSITORY https://github.com/CLIUtils/CLI11.git
GIT_TAG v2.5.0
)
FetchContent_MakeAvailable(cli11)
# add_library(poetlib
# Base/Grid.cpp
# Base/SimParams.cpp
@ -86,16 +70,16 @@ target_compile_definitions(POETLib PUBLIC STRICT_R_HEADERS OMPI_SKIP_MPICXX)
file(READ "${PROJECT_SOURCE_DIR}/R_lib/kin_r_library.R" R_KIN_LIB )
file(READ "${PROJECT_SOURCE_DIR}/R_lib/init_r_lib.R" R_INIT_LIB)
file(READ "${PROJECT_SOURCE_DIR}/R_lib/ai_surrogate_model.R" R_AI_SURROGATE_LIB)
configure_file(poet.hpp.in poet.hpp @ONLY)
add_executable(poet poet.cpp)
target_link_libraries(poet PRIVATE POETLib MPI::MPI_C RRuntime CLI11::CLI11)
target_link_libraries(poet PRIVATE POETLib MPI::MPI_C RRuntime)
target_include_directories(poet PRIVATE "${CMAKE_CURRENT_BINARY_DIR}")
add_executable(poet_init initializer.cpp)
target_link_libraries(poet_init PRIVATE POETLib RRuntime CLI11::CLI11)
target_link_libraries(poet_init PRIVATE POETLib RRuntime)
target_include_directories(poet_init PRIVATE "${CMAKE_CURRENT_BINARY_DIR}")
install(TARGETS poet poet_init DESTINATION bin)

9
src/Chemistry/ChemistryDefs.hpp
View File

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

97
src/Chemistry/ChemistryModule.cpp
View File

@ -1,14 +1,11 @@
#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>
@ -66,8 +63,7 @@ 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] = distance != 0 ? 1 / std::sqrt(distance) : 0;
weights[point_i] = 1 / std::sqrt(distance);
assert(!std::isnan(weights[point_i]));
inv_sum += weights[point_i];
}
@ -98,9 +94,63 @@ 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;
}
@ -117,12 +167,16 @@ poet::ChemistryModule::ChemistryModule(
this->n_cells = chem_params.total_grid_cells;
if (!is_master) {
PhreeqcMatrix pqc_mat =
PhreeqcMatrix(chem_params.database, chem_params.pqc_script,
chem_params.with_h0_o0, chem_params.with_redox);
this->pqc_runner =
std::make_unique<PhreeqcRunner>(pqc_mat.subset(chem_params.pqc_ids));
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_);
// }
}
}
@ -133,10 +187,11 @@ poet::ChemistryModule::~ChemistryModule() {
}
void poet::ChemistryModule::initializeDHT(
uint32_t size_mb, const NamedVector<std::uint32_t> &key_species,
bool has_het_ids) {
uint32_t size_mb, const NamedVector<std::uint32_t> &key_species) {
constexpr uint32_t MB_FACTOR = 1E6;
this->dht_enabled = true;
MPI_Comm dht_comm;
if (this->is_master) {
@ -166,7 +221,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, has_het_ids);
this->prop_count, interp_enabled);
this->dht->setBaseTotals(base_totals.at(0), base_totals.at(1));
}
}
@ -257,10 +312,9 @@ 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;
}
}
@ -317,8 +371,3 @@ void poet::ChemistryModule::unshuffleField(const std::vector<double> &in_buffer,
}
}
}
void poet::ChemistryModule::set_ai_surrogate_validity_vector(
std::vector<int> r_vector) {
this->ai_surrogate_validity_vector = r_vector;
}

35
src/Chemistry/ChemistryModule.hpp
View File

@ -8,11 +8,10 @@
#include "ChemistryDefs.hpp"
#include "Init/InitialList.hpp"
#include "NameDouble.h"
#include "SurrogateModels/DHT_Wrapper.hpp"
#include "SurrogateModels/Interpolation.hpp"
#include "PhreeqcRunner.hpp"
#include "PhreeqcEngine.hpp"
#include <array>
#include <cstdint>
#include <map>
@ -76,19 +75,14 @@ public:
struct SurrogateSetup {
std::vector<std::string> prop_names;
std::array<double, 2> base_totals;
bool has_het_ids;
bool dht_enabled;
std::uint32_t dht_size_mb;
int dht_snaps;
std::string dht_out_dir;
bool interp_enabled;
std::uint32_t interp_bucket_size;
std::uint32_t interp_size_mb;
std::uint32_t interp_min_entries;
bool ai_surrogate_enabled;
};
void masterEnableSurrogates(const SurrogateSetup &setup) {
@ -98,17 +92,9 @@ public:
this->dht_enabled = setup.dht_enabled;
this->interp_enabled = setup.interp_enabled;
this->ai_surrogate_enabled = setup.ai_surrogate_enabled;
this->base_totals = setup.base_totals;
if (this->dht_enabled || this->interp_enabled) {
this->initializeDHT(setup.dht_size_mb, this->params.dht_species,
setup.has_het_ids);
if (setup.dht_snaps != DHT_SNAPS_DISABLED) {
this->setDHTSnapshots(setup.dht_snaps, setup.dht_out_dir);
}
this->initializeDHT(setup.dht_size_mb, this->params.dht_species);
}
if (this->interp_enabled) {
@ -233,11 +219,6 @@ public:
this->print_progessbar = enabled;
};
/**
* **Master only** Set the ai surrogate validity vector from R
*/
void set_ai_surrogate_validity_vector(std::vector<int> r_vector);
std::vector<uint32_t> GetWorkerInterpolationCalls() const;
std::vector<double> GetWorkerInterpolationWriteTimings() const;
@ -247,12 +228,9 @@ public:
std::vector<uint32_t> GetWorkerPHTCacheHits() const;
std::vector<int> ai_surrogate_validity_vector;
protected:
void initializeDHT(uint32_t size_mb,
const NamedVector<std::uint32_t> &key_species,
bool has_het_ids);
const NamedVector<std::uint32_t> &key_species);
void setDHTSnapshots(int type, const std::string &out_dir);
void setDHTReadFile(const std::string &input_file);
@ -271,8 +249,7 @@ protected:
CHEM_IP_SIGNIF_VEC,
CHEM_WORK_LOOP,
CHEM_PERF,
CHEM_BREAK_MAIN_LOOP,
CHEM_AI_BCAST_VALIDITY
CHEM_BREAK_MAIN_LOOP
};
enum { LOOP_WORK, LOOP_END };
@ -371,8 +348,6 @@ protected:
bool interp_enabled{false};
std::unique_ptr<poet::InterpolationModule> interp;
bool ai_surrogate_enabled{false};
static constexpr uint32_t BUFFER_OFFSET = 5;
inline void ChemBCast(void *buf, int count, MPI_Datatype datatype) const {
@ -403,7 +378,7 @@ protected:
const InitialList::ChemistryInit params;
std::unique_ptr<PhreeqcRunner> pqc_runner;
std::map<int, std::unique_ptr<PhreeqcEngine>> phreeqc_instances;
};
} // namespace poet

38
src/Chemistry/MasterFunctions.cpp
View File

@ -159,20 +159,6 @@ std::vector<uint32_t> poet::ChemistryModule::GetWorkerPHTCacheHits() const {
return ret;
}
inline std::vector<int> shuffleVector(const std::vector<int> &in_vector,
uint32_t size_per_prop,
uint32_t wp_count) {
std::vector<int> out_buffer(in_vector.size());
uint32_t write_i = 0;
for (uint32_t i = 0; i < wp_count; i++) {
for (uint32_t j = i; j < size_per_prop; j += wp_count) {
out_buffer[write_i] = in_vector[j];
write_i++;
}
}
return out_buffer;
}
inline std::vector<double> shuffleField(const std::vector<double> &in_field,
uint32_t size_per_prop,
uint32_t prop_count,
@ -261,10 +247,8 @@ inline void poet::ChemistryModule::MasterSendPkgs(
send_buffer[end_of_wp + 2] = dt;
// current time of simulation (age) in seconds
send_buffer[end_of_wp + 3] = this->simtime;
// current work package start location in field
uint32_t wp_start_index = std::accumulate(wp_sizes_vector.begin(), std::next(wp_sizes_vector.begin(), count_pkgs), 0);
send_buffer[end_of_wp + 4] = wp_start_index;
// placeholder for work_package_count
send_buffer[end_of_wp + 4] = 0.;
/* ATTENTION Worker p has rank p+1 */
// MPI_Send(send_buffer, end_of_wp + BUFFER_OFFSET, MPI_DOUBLE, p + 1,
@ -368,21 +352,8 @@ void poet::ChemistryModule::MasterRunParallel(double dt) {
int pkg_to_send, pkg_to_recv;
int free_workers;
int i_pkgs;
int ftype;
const std::vector<uint32_t> wp_sizes_vector =
CalculateWPSizesVector(this->n_cells, this->wp_size);
if (this->ai_surrogate_enabled) {
ftype = CHEM_AI_BCAST_VALIDITY;
PropagateFunctionType(ftype);
this->ai_surrogate_validity_vector = shuffleVector(this->ai_surrogate_validity_vector,
this->n_cells,
wp_sizes_vector.size());
ChemBCast(&this->ai_surrogate_validity_vector.front(), this->n_cells, MPI_INT);
}
ftype = CHEM_WORK_LOOP;
int ftype = CHEM_WORK_LOOP;
PropagateFunctionType(ftype);
MPI_Barrier(this->group_comm);
@ -392,6 +363,9 @@ void poet::ChemistryModule::MasterRunParallel(double dt) {
/* start time measurement of sequential part */
seq_a = MPI_Wtime();
const std::vector<uint32_t> wp_sizes_vector =
CalculateWPSizesVector(this->n_cells, this->wp_size);
/* shuffle grid */
// grid.shuffleAndExport(mpi_buffer);
std::vector<double> mpi_buffer =

44
src/Chemistry/SurrogateModels/DHT.c
View File

@ -1,3 +1,4 @@
/// Time-stamp: "Last modified 2023-08-10 11:50:46 mluebke"
/*
** Copyright (C) 2017-2021 Max Luebke (University of Potsdam)
**
@ -105,8 +106,6 @@ 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
@ -191,9 +190,6 @@ 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);
@ -284,9 +280,6 @@ 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],
@ -558,41 +551,6 @@ 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;

48
src/Chemistry/SurrogateModels/DHT.h
View File

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

89
src/Chemistry/SurrogateModels/DHT_Wrapper.cpp
View File

@ -1,3 +1,4 @@
// Time-stamp: "Last modified 2023-11-01 10:54:45 mluebke"
/*
** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
@ -24,7 +25,6 @@
#include "Init/InitialList.hpp"
#include "Rounding.hpp"
#include <Rcpp/proxy/ProtectedProxy.h>
#include <algorithm>
#include <cassert>
#include <cmath>
@ -43,12 +43,11 @@ 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,
bool _has_het_ids)
uint32_t data_count, bool _with_interp)
: 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), has_het_ids(_has_het_ids) {
with_interp(_with_interp) {
// initialize DHT object
// key size = count of key elements + timestep
uint32_t key_size = (key_count + 1) * sizeof(Lookup_Keyelement);
@ -129,47 +128,42 @@ 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) {
continue;
}
if (work_package.mapping[i] == CHEM_PQC) {
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]);
// 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;
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;
}
uint32_t proc, index;
auto &key = dht_results.keys[i];
const auto data =
(with_interp ? outputToInputAndRates(work_package.input[i],
work_package.output[i])
: work_package.output[i]);
// void *data = (void *)&(work_package[i * this->data_count]);
// fuzz data (round, logarithm etc.)
// insert simulated data with fuzzed key into DHT
int res = DHT_write(this->dht_object, key.data(),
const_cast<double *>(data.data()), &proc, &index);
dht_results.locations[i] = {proc, index};
// if data was successfully written ...
if ((res != DHT_SUCCESS) && (res == DHT_WRITE_SUCCESS_WITH_EVICTION)) {
dht_evictions++;
}
dht_results.filledDHT[i] = true;
}
}
@ -273,10 +267,9 @@ LookupKey DHT_Wrapper::fuzzForDHT_R(const std::vector<double> &cell,
NamedVector<double> input_nv(this->output_names, cell);
const std::vector<double> eval_vec =
Rcpp::as<std::vector<double>>(hooks.dht_fuzz(input_nv));
const std::vector<double> eval_vec = hooks.dht_fuzz(input_nv);
assert(eval_vec.size() == this->key_count);
LookupKey vecFuzz(this->key_count + 1 + has_het_ids, {.0});
LookupKey vecFuzz(this->key_count + 1, {.0});
DHT_Rounder rounder;
@ -296,9 +289,6 @@ 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;
}
@ -306,7 +296,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 + has_het_ids, {.0});
LookupKey vecFuzz(this->key_count + 1, {.0});
DHT_Rounder rounder;
int totals_i = 0;
@ -332,9 +322,6 @@ 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;
}

4
src/Chemistry/SurrogateModels/DHT_Wrapper.hpp
View File

@ -1,3 +1,4 @@
// Time-stamp: "Last modified 2023-09-08 14:43:02 mluebke"
/*
** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
@ -87,7 +88,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, bool has_het_ids);
uint32_t data_count, bool with_interp);
/**
* @brief Destroy the dht wrapper object
*
@ -264,7 +265,6 @@ private:
DHT_ResultObject dht_results;
std::array<double, 2> base_totals{0};
bool has_het_ids{false};
};
} // namespace poet

1
src/Chemistry/SurrogateModels/HashFunctions.cpp
View File

@ -1,3 +1,4 @@
// Time-stamp: "Last modified 2023-04-24 23:20:55 mluebke"
/*
**-----------------------------------------------------------------------------
** MurmurHash2 was written by Austin Appleby, and is placed in the public

1
src/Chemistry/SurrogateModels/HashFunctions.hpp
View File

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

5
src/Chemistry/SurrogateModels/Interpolation.hpp
View File

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

30
src/Chemistry/SurrogateModels/InterpolationModule.cpp
View File

@ -1,3 +1,4 @@
// Time-stamp: "Last modified 2023-08-16 17:02:31 mluebke"
#include "Init/InitialList.hpp"
#include "Interpolation.hpp"
@ -8,13 +9,11 @@
#include "Rounding.hpp"
#include <Rcpp.h>
#include <Rcpp/proxy/ProtectedProxy.h>
#include <Rinternals.h>
#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <functional>
@ -76,11 +75,6 @@ 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;
@ -100,8 +94,7 @@ void InterpolationModule::tryInterpolation(WorkPackage &work_package) {
if (hooks.interp_pre.isValid()) {
NamedVector<double> nv_in(this->out_names, work_package.input[wp_i]);
std::vector<int> rm_indices = Rcpp::as<std::vector<int>>(
hooks.interp_pre(nv_in, pht_result.in_values));
auto rm_indices = hooks.interp_pre(nv_in, pht_result.in_values);
pht_result.size -= rm_indices.size();
@ -122,30 +115,15 @@ 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(to_calc_cache[cell_id], work_package.input[wp_i],
f_interpolate(dht_instance.getKeyElements(), 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 (Rcpp::as<bool>(hooks.interp_post(nv_result))) {
if (hooks.interp_post(nv_result)) {
interp_result.status[wp_i] = INSUFFICIENT_DATA;
continue;
}

17
src/Chemistry/SurrogateModels/LookupKey.hpp
View File

@ -1,3 +1,4 @@
// Time-stamp: "Last modified 2023-08-11 10:12:52 mluebke"
#ifndef LOOKUPKEY_H_
#define LOOKUPKEY_H_
@ -10,21 +11,9 @@
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 {
@ -35,10 +24,6 @@ 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> {

1
src/Chemistry/SurrogateModels/ProximityHashTable.cpp
View File

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

13
src/Chemistry/SurrogateModels/Rounding.hpp
View File

@ -20,11 +20,6 @@ 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))));
@ -65,14 +60,6 @@ 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)))) -

64
src/Chemistry/WorkerFunctions.cpp
View File

@ -9,6 +9,7 @@
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <map>
#include <mpi.h>
#include <stdexcept>
#include <string>
@ -46,16 +47,6 @@ void poet::ChemistryModule::WorkerLoop() {
switch (func_type) {
case CHEM_FIELD_INIT: {
ChemBCast(&this->prop_count, 1, MPI_UINT32_T);
if (this->ai_surrogate_enabled) {
this->ai_surrogate_validity_vector.resize(
this->n_cells); // resize statt reserve?
}
break;
}
case CHEM_AI_BCAST_VALIDITY: {
// Receive the index vector of valid ai surrogate predictions
MPI_Bcast(&this->ai_surrogate_validity_vector.front(), this->n_cells,
MPI_INT, 0, this->group_comm);
break;
}
case CHEM_WORK_LOOP: {
@ -127,7 +118,7 @@ void poet::ChemistryModule::WorkerDoWork(MPI_Status &probe_status,
uint32_t iteration;
double dt;
double current_sim_time;
uint32_t wp_start_index;
int count = double_count;
std::vector<double> mpi_buffer(count);
@ -153,8 +144,8 @@ void poet::ChemistryModule::WorkerDoWork(MPI_Status &probe_status,
// current simulation time ('age' of simulation)
current_sim_time = mpi_buffer[count + 3];
// current work package start location in field
wp_start_index = mpi_buffer[count + 4];
/* 4th double value is currently a placeholder */
// placeholder = mpi_buffer[count+4];
for (std::size_t wp_i = 0; wp_i < s_curr_wp.size; wp_i++) {
s_curr_wp.input[wp_i] =
@ -179,15 +170,6 @@ void poet::ChemistryModule::WorkerDoWork(MPI_Status &probe_status,
interp->tryInterpolation(s_curr_wp);
}
if (this->ai_surrogate_enabled) {
// Map valid predictions from the ai surrogate in the workpackage
for (int i = 0; i < s_curr_wp.size; i++) {
if (this->ai_surrogate_validity_vector[wp_start_index + i] == 1) {
s_curr_wp.mapping[i] = CHEM_AISURR;
}
}
}
phreeqc_time_start = MPI_Wtime();
WorkerRunWorkPackage(s_curr_wp, current_sim_time, dt);
@ -246,17 +228,6 @@ 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()");
@ -311,19 +282,28 @@ 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) {
to_ignore.push_back(wp_id);
continue;
}
}
this->pqc_runner->run(inout_chem, dTimestep, to_ignore);
for (std::size_t wp_id = 0; wp_id < work_package.size; wp_id++) {
if (work_package.mapping[wp_id] == CHEM_PQC) {
work_package.output[wp_id] = inout_chem[wp_id];
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++];
}
}
}
}

52
src/Init/ChemistryInit.cpp
View File

@ -2,14 +2,16 @@
#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"]);
}
@ -30,30 +32,6 @@ void InitialList::initChemistry(const Rcpp::List &chem) {
}
}
if (chem.containsElementNamed("ai_surrogate_input_script")) {
std::string ai_surrogate_input_script_path = chem["ai_surrogate_input_script"];
ai_surrogate_input_script_path = Rcpp::as<std::string>(Rcpp::Function("normalizePath")(Rcpp::wrap(ai_surrogate_input_script_path)));
// Copying the entire script for the init file
std::ifstream file(ai_surrogate_input_script_path);
if (!file.is_open()) {
// print error message and return
Rcpp::Rcerr << "AI surrogate input script was not found at: " << ai_surrogate_input_script_path << std::endl;
}
std::stringstream buffer;
buffer << file.rdbuf();
std::string fileContent = buffer.str();
file.close();
// Get base path
ai_surrogate_input_script_path = ai_surrogate_input_script_path.substr(0, ai_surrogate_input_script_path.find_last_of('/') + 1);
// Add the filepath as a global variable in R to enable relative filepaths in the R script
fileContent += "\nai_surrogate_base_path <- \"" + ai_surrogate_input_script_path + "\"";
this->ai_surrogate_input_script = fileContent;
}
this->field_header =
Rcpp::as<std::vector<std::string>>(this->initial_grid.names());
this->field_header.erase(this->field_header.begin());
@ -67,22 +45,26 @@ 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_ids.size(); ++i) {
// chem_init.pqc_input[pqc_ids[i]] = pqc_scripts[i];
// }
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};
}
// chem_init.pqc_sol_order = pqc_solutions;
chem_init.dht_species = dht_species;
chem_init.interp_species = interp_species;
chem_init.ai_surrogate_input_script = ai_surrogate_input_script;
if (this->chem_hooks.size() > 0) {
if (this->chem_hooks.containsElementNamed("dht_fill")) {

178
src/Init/DiffusionInit.cpp
View File

@ -1,3 +1,4 @@
#include <algorithm>
#include <tug/Boundary.hpp>
// leave above Rcpp includes, as eigen seem to have problems with a preceding
// Rcpp include
@ -7,7 +8,7 @@
#include "DataStructures/Field.hpp"
#include "InitialList.hpp"
#include "PhreeqcMatrix.hpp"
#include "PhreeqcInit.hpp"
#include <Rcpp/Function.h>
#include <Rcpp/proxy/ProtectedProxy.h>
@ -26,9 +27,9 @@ namespace poet {
enum SEXP_TYPE { SEXP_IS_LIST, SEXP_IS_VEC };
const std::map<std::uint8_t, std::string> tug_side_mapping = {
{tug::BC_SIDE_RIGHT , "E"},
{tug::BC_SIDE_LEFT , "W"},
{tug::BC_SIDE_TOP , "N"},
{tug::BC_SIDE_RIGHT, "E"},
{tug::BC_SIDE_LEFT, "W"},
{tug::BC_SIDE_TOP, "N"},
{tug::BC_SIDE_BOTTOM, "S"}};
static std::vector<TugType> colMajToRowMaj(const Rcpp::NumericVector &vec,
@ -52,104 +53,100 @@ 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 PhreeqcMatrix &phreeqc_mat) {
const Rcpp::List &inner_boundaries) {
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) {
for (const auto &species : this->transport_names) {
Rcpp::List spec_list;
for (const auto &side : tug_side_mapping) {
@ -166,8 +163,8 @@ InitialList::resolveBoundaries(const Rcpp::List &boundaries_list,
if (boundaries_list.containsElementNamed(side.second.c_str())) {
const Rcpp::List mapping = boundaries_list[side.second];
const Rcpp::IntegerVector cells = mapping["cell"]; // MDL 2024-06-13
const Rcpp::IntegerVector values = mapping["sol_id"]; // MDL
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()) {
@ -182,11 +179,8 @@ 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] =
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]));
c_value[c_id] = Rcpp::as<TugType>(
resolve_R(this->phreeqc_mat, Rcpp::wrap(species), values[i]));
} else {
throw std::runtime_error("Unknown boundary type");
}
@ -197,7 +191,7 @@ InitialList::resolveBoundaries(const Rcpp::List &boundaries_list,
Rcpp::Named("type") = c_type, Rcpp::Named("value") = c_value);
}
bound_list[species_name] = spec_list;
bound_list[species] = spec_list;
if (inner_boundaries.size() > 0) {
const Rcpp::NumericVector &inner_row_vec = inner_boundaries["row"];
@ -217,14 +211,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(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])));
c_value.push_back(Rcpp::as<TugType>(resolve_R(
this->phreeqc_mat, Rcpp::wrap(species), inner_pqc_id_vec[i])));
}
inner_bound[species_name] =
inner_bound[species] =
Rcpp::List::create(Rcpp::Named("row") = Rcpp::wrap(rows),
Rcpp::Named("col") = Rcpp::wrap(cols),
Rcpp::Named("value") = Rcpp::wrap(c_value));
@ -285,8 +276,7 @@ static Rcpp::List parseAlphas(const SEXP &input,
return out_list;
}
void InitialList::initDiffusion(const Rcpp::List &diffusion_input,
const PhreeqcMatrix &phreeqc) {
void InitialList::initDiffusion(const Rcpp::List &diffusion_input) {
Rcpp::List boundaries;
Rcpp::List inner_boundaries;
@ -308,7 +298,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, phreeqc);
resolveBoundaries(boundaries, inner_boundaries);
this->boundaries = resolved_boundaries.first;
this->inner_boundaries = resolved_boundaries.second;
@ -413,4 +403,4 @@ InitialList::DiffusionInit InitialList::getDiffusionInit() const {
return diff_init;
}
} // namespace poet
} // namespace poet

184
src/Init/GridInit.cpp
View File

@ -1,19 +1,13 @@
// 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 "PhreeqcMatrix.hpp"
#include "PhreeqcInit.hpp"
#include <RInside.h>
#include <Rcpp/Function.h>
#include <Rcpp/vector/Matrix.h>
#include <Rcpp/vector/instantiation.h>
#include <fstream>
#include <cstdint>
#include <map>
#include <memory>
#include <regex>
#include <sstream>
#include <string>
@ -21,6 +15,38 @@
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) {
@ -34,6 +60,25 @@ 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');
@ -53,26 +98,10 @@ static std::string readFile(const std::string &path) {
return buffer.str();
}
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) {
void InitialList::prepareGrid(const Rcpp::List &grid_input) {
// parse input values
Rcpp::Function unique_R("unique");
std::string script;
std::string database;
@ -97,9 +126,8 @@ PhreeqcMatrix InitialList::prepareGrid(const Rcpp::List &grid_input) {
}
this->database = database;
this->pqc_script = script;
Rcpp::IntegerMatrix grid_def =
Rcpp::NumericMatrix grid_def =
grid_input[GRID_MEMBER_STR(GridMembers::GRID_DEF)];
Rcpp::NumericVector grid_size =
grid_input[GRID_MEMBER_STR(GridMembers::GRID_SIZE)];
@ -128,77 +156,35 @@ PhreeqcMatrix InitialList::prepareGrid(const Rcpp::List &grid_input) {
throw std::runtime_error("Grid size must be positive.");
}
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 = std::make_unique<PhreeqcInit>(database, script);
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;
this->phreeqc_mat = pqcScriptToGrid(phreeqc, R);
this->initial_grid = matToGrid(R, this->phreeqc_mat, grid_def);
if (with_h0_o0 && !with_redox) {
throw std::runtime_error(
"Output of H(0) and O(0) can only be used with redox.");
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));
}
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

80
src/Init/InitialList.cpp
View File

@ -1,6 +1,5 @@
#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>
@ -11,8 +10,8 @@
namespace poet {
void InitialList::initializeFromList(const Rcpp::List &setup) {
PhreeqcMatrix phreeqc = prepareGrid(setup[grid_key]);
initDiffusion(setup[diffusion_key], phreeqc);
prepareGrid(setup[grid_key]);
initDiffusion(setup[diffusion_key]);
initChemistry(setup[chemistry_key]);
}
@ -54,30 +53,26 @@ 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)]);
@ -87,8 +82,6 @@ void InitialList::importList(const Rcpp::List &setup, bool minimal) {
this->chem_hooks =
Rcpp::as<Rcpp::List>(setup[static_cast<int>(ExportList::CHEM_HOOKS)]);
this->ai_surrogate_input_script = Rcpp::as<std::string>(setup[static_cast<int>(ExportList::AI_SURROGATE_INPUT_SCRIPT)]);
}
Rcpp::List InitialList::exportList() {
@ -113,34 +106,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);
out[static_cast<int>(ExportList::CHEM_HOOKS)] = this->chem_hooks;
out[static_cast<int>(ExportList::AI_SURROGATE_INPUT_SCRIPT)] = this->ai_surrogate_input_script;
return out;
}

76
src/Init/InitialList.hpp
View File

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

109
src/Transport/DiffusionModule.cpp
View File

@ -1,48 +1,97 @@
// SPDX-License-Identifier: GPL-2.0-or-later
/*
* DiffusionModule.cpp - Coupling module betweeen POET and tug
* Copyright (C) 2018-2025 Max Luebke (University of Potsdam)
*/
** 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.
*/
#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 <tug/Boundary.hpp>
#include <tug/Diffusion.hpp>
#include <vector>
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) {
#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
auto &species_conc = this->transport_field[sol_name];
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();
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);
boundary.deserialize(this->param_list.boundaries[sol_name]);
@ -50,8 +99,14 @@ void DiffusionModule::simulate(double requested_dt) {
boundary.setInnerBoundaries(this->param_list.inner_boundaries[sol_name]);
}
diffusion_solver.setTimestep(requested_dt);
diffusion_solver.run();
grid.setAlpha(alpha_x, alpha_y);
grid.setConcentrations(conc);
sim.setTimestep(requested_dt);
sim.run();
species_conc = MatrixToVec(grid.getConcentrations());
}
const auto end_diffusion_t = std::chrono::high_resolution_clock::now();

81
src/initializer.cpp
View File

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

443
src/poet.cpp
View File

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

67
src/poet.hpp.in
View File

@ -23,6 +23,7 @@
#pragma once
#include <cstdint>
#include <set>
#include <string>
#include <vector>
@ -34,39 +35,61 @@ static const char *poet_version = "@POET_VERSION@";
static const inline std::string kin_r_library = R"(@R_KIN_LIB@)";
static const inline std::string init_r_library = R"(@R_INIT_LIB@)";
static const inline std::string ai_surrogate_r_library =
R"(@R_AI_SURROGATE_LIB@)";
static const inline std::string r_runtime_parameters = "mysetup";
const std::set<std::string> flaglist{"ignore-result", "dht", "P", "progress",
"interp"};
const std::set<std::string> paramlist{
"work-package-size", "dht-strategy", "dht-size", "dht-snaps",
"dht-file", "interp-size", "interp-min", "interp-bucket-entries"};
constexpr uint32_t CHEM_DEFAULT_WORK_PACKAGE_SIZE = 32;
constexpr uint32_t CHEM_DHT_SIZE_PER_PROCESS_MB = 1.5E3;
struct RuntimeParameters {
std::string out_dir;
std::vector<double> timesteps;
bool print_progressbar;
uint32_t work_package_size;
Rcpp::List init_params;
// MDL added to accomodate for qs::qsave/qread
bool as_rds = false;
bool as_qs = false;
std::string out_ext;
bool use_dht;
std::uint32_t dht_size;
std::uint8_t dht_snaps;
bool print_progress = false;
bool use_interp;
std::uint32_t interp_size;
std::uint32_t interp_min_entries;
std::uint32_t interp_bucket_entries;
static constexpr std::uint32_t WORK_PACKAGE_SIZE_DEFAULT = 32;
std::uint32_t work_package_size = WORK_PACKAGE_SIZE_DEFAULT;
struct ChemistryParams {
// std::string database_path;
// std::string input_script;
bool use_dht = false;
static constexpr std::uint32_t DHT_SIZE_DEFAULT = 1.5E3;
std::uint32_t dht_size = DHT_SIZE_DEFAULT;
static constexpr std::uint8_t DHT_SNAPS_DEFAULT = 0;
std::uint8_t dht_snaps = DHT_SNAPS_DEFAULT;
// bool use_dht;
// std::uint64_t dht_size;
// int dht_snaps;
// std::string dht_file;
// std::string dht_outdir;
// NamedVector<std::uint32_t> dht_signifs;
bool use_interp = false;
static constexpr std::uint32_t INTERP_SIZE_DEFAULT = 100;
std::uint32_t interp_size = INTERP_SIZE_DEFAULT;
static constexpr std::uint32_t INTERP_MIN_ENTRIES_DEFAULT = 5;
std::uint32_t interp_min_entries = INTERP_MIN_ENTRIES_DEFAULT;
static constexpr std::uint32_t INTERP_BUCKET_ENTRIES_DEFAULT = 20;
std::uint32_t interp_bucket_entries = INTERP_BUCKET_ENTRIES_DEFAULT;
// bool use_interp;
// std::uint64_t pht_size;
// std::uint32_t pht_max_entries;
// std::uint32_t interp_min_entries;
// NamedVector<std::uint32_t> pht_signifs;
bool use_ai_surrogate = false;
// struct Chem_Hook_Functions {
// RHookFunction<bool> dht_fill;
// RHookFunction<std::vector<double>> dht_fuzz;
// RHookFunction<std::vector<std::size_t>> interp_pre;
// RHookFunction<bool> interp_post;
// } hooks;
// void initFromR(RInsidePOET &R);
};
};

4
test/testField.cpp
View File

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

10
test/testNamedVector.cpp
View File

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