mirror of
https://git.gfz-potsdam.de/naaice/iphreeqc.git
synced 2025-12-16 16:44:49 +01:00
Merge branch 'ml/improve-API' into 'poet'
Improve usability outside POET See merge request naaice/iphreeqc!10
This commit is contained in:
Marco De Lucia
commit
9061e210e6
148
.gitlab-ci.yml
148
.gitlab-ci.yml
@ -1,134 +1,22 @@
|
|||||||
#
|
image: ubuntu:24.04
|
||||||
# https://code.chs.usgs.gov/coupled/iphreeqc
|
|
||||||
# SRC 2020-12-02T18:39:55-07:00
|
|
||||||
# SRC 2021-10-31T11:23:46-06:00 -- changed pull to squash -- HEAD:9499e78cb2493da6f56683ea669cd7ed23541ddc
|
|
||||||
#
|
|
||||||
image: ${CI_REGISTRY}/coupled/containers/buildpack-deps:bionic-scm
|
|
||||||
|
|
||||||
stages:
|
|
||||||
- sync
|
|
||||||
- trigger
|
|
||||||
|
|
||||||
before_script:
|
before_script:
|
||||||
- eval $(ssh-agent -s)
|
- apt-get update -y
|
||||||
- echo "${SSH_PRIVATE_KEY_ENC}" | base64 --decode | tr -d '\r' | ssh-add -
|
- apt-get install -y build-essential cmake git
|
||||||
- mkdir -p ~/.ssh
|
|
||||||
- chmod 700 ~/.ssh
|
|
||||||
- ssh-keyscan ${CI_SERVER_HOST} >> ~/.ssh/known_hosts
|
|
||||||
- chmod 644 ~/.ssh/known_hosts
|
|
||||||
- git config --global user.email "darth@empire.com"
|
|
||||||
- git config --global user.name "Darth Vader"
|
|
||||||
|
|
||||||
subtree-sync:
|
stages:
|
||||||
stage: sync
|
- test
|
||||||
|
|
||||||
##
|
test:
|
||||||
## Only run if on the master branch and the variable GROUP is set
|
stage: test
|
||||||
##
|
script:
|
||||||
## change this to
|
- mkdir _build && cd _build
|
||||||
## only:
|
- cmake -DBUILD_SHARED_LIBS=ON -DCMAKE_BUILD_TYPE=Release ..
|
||||||
## - master@$GROUP/iphreeqc
|
- make -j$(nproc)
|
||||||
## and set GROUP to coupled before merge
|
- ctest --output-junit test_results.xml
|
||||||
only:
|
artifacts:
|
||||||
refs:
|
when: always
|
||||||
- master
|
paths:
|
||||||
variables:
|
- _build/test_results.xml
|
||||||
- $GROUP
|
reports:
|
||||||
|
junit: _build/test_results.xml
|
||||||
script:
|
|
||||||
##
|
|
||||||
## Must re-clone in order for the subtree merge to work
|
|
||||||
## tried re-setting the url for the origin but didn't work
|
|
||||||
##
|
|
||||||
- cd ..
|
|
||||||
- rm -rf ${CI_PROJECT_NAME}
|
|
||||||
- git clone git@${CI_SERVER_HOST}:${CI_PROJECT_PATH}.git
|
|
||||||
- cd ${CI_PROJECT_NAME}
|
|
||||||
|
|
||||||
##
|
|
||||||
## Sync subtrees
|
|
||||||
##
|
|
||||||
- |
|
|
||||||
#!/bin/bash -ex
|
|
||||||
#
|
|
||||||
# iphreeqc/ git@${CI_SERVER_HOST}:${GROUP}/iphreeqc.git
|
|
||||||
# ├─database/ ├─git@${CI_SERVER_HOST}:${GROUP}/subtrees/phreeqc3-database.git* database*
|
|
||||||
# ├─examples/ │ └─examples
|
|
||||||
# │ ├─c/ │ ├─git@${CI_SERVER_HOST}:${GROUP}/subtrees/phreeqc-commanuscript-cgfinal-examples-c.git* examples/c*
|
|
||||||
# │ ├─com/ │ ├─git@${CI_SERVER_HOST}:${GROUP}/subtrees/phreeqc-commanuscript-cgfinal-examples-com.git* examples/com*
|
|
||||||
# │ └─fortran/ │ └─git@${CI_SERVER_HOST}:${GROUP}/subtrees/phreeqc-COMManuscript-CGfinal-examples-fortran.git* examples/fortran*
|
|
||||||
# ├─phreeqc3-doc/ ├─git@${CI_SERVER_HOST}:${GROUP}/subtrees/phreeqc3-doc.git* phreeqc3-doc*
|
|
||||||
# ├─phreeqc3-examples/ ├─git@${CI_SERVER_HOST}:${GROUP}/subtrees/phreeqc3-examples.git* phreeqc3-examples*
|
|
||||||
# └─src/ └─git@${CI_SERVER_HOST}:${GROUP}/subtrees/iphreeqc-src.git% src%
|
|
||||||
# └─phreeqcpp/ └─git@${CI_SERVER_HOST}:${GROUP}/subtrees/phreeqc3-src.git src/phreeqcpp
|
|
||||||
# └─common/ └─git@${CI_SERVER_HOST}:${GROUP}/subtrees/phreeqc3-src-common.git src/phreeqcpp/common
|
|
||||||
|
|
||||||
git_subtree() {
|
|
||||||
git subtree "${1}" --prefix="${2}" "${4}" master 2>&1 | grep -v "^[[:digit:]].*/[[:digit:]].*"
|
|
||||||
}
|
|
||||||
|
|
||||||
declare -A urls=( \
|
|
||||||
["iphreeqc-src"]="git@${CI_SERVER_HOST}:${GROUP}/subtrees/iphreeqc-src.git" \
|
|
||||||
["phreeqc-commanuscript-cgfinal-examples-c"]="git@${CI_SERVER_HOST}:${GROUP}/subtrees/phreeqc-commanuscript-cgfinal-examples-c.git" \
|
|
||||||
["phreeqc-commanuscript-cgfinal-examples-com"]="git@${CI_SERVER_HOST}:${GROUP}/subtrees/phreeqc-commanuscript-cgfinal-examples-com.git" \
|
|
||||||
["phreeqc-COMManuscript-CGfinal-examples-fortran"]="git@${CI_SERVER_HOST}:${GROUP}/subtrees/phreeqc-COMManuscript-CGfinal-examples-fortran.git" \
|
|
||||||
["phreeqc3-database"]="git@${CI_SERVER_HOST}:${GROUP}/subtrees/phreeqc3-database.git" \
|
|
||||||
["phreeqc3-doc"]="git@${CI_SERVER_HOST}:${GROUP}/subtrees/phreeqc3-doc.git" \
|
|
||||||
["phreeqc3-examples"]="git@${CI_SERVER_HOST}:${GROUP}/subtrees/phreeqc3-examples.git" \
|
|
||||||
)
|
|
||||||
|
|
||||||
declare -A prefixes=( \
|
|
||||||
["iphreeqc-src"]="src" \
|
|
||||||
["phreeqc-commanuscript-cgfinal-examples-c"]="examples/c" \
|
|
||||||
["phreeqc-commanuscript-cgfinal-examples-com"]="examples/com" \
|
|
||||||
["phreeqc-COMManuscript-CGfinal-examples-fortran"]="examples/fortran" \
|
|
||||||
["phreeqc3-database"]="database" \
|
|
||||||
["phreeqc3-doc"]="phreeqc3-doc" \
|
|
||||||
["phreeqc3-examples"]="phreeqc3-examples" \
|
|
||||||
)
|
|
||||||
|
|
||||||
export GIT_EDITOR=true
|
|
||||||
|
|
||||||
for remote in "${!urls[@]}"; do
|
|
||||||
# git_subtree "pull" "${prefixes[$remote]}" "$remote" "${urls[$remote]}"
|
|
||||||
git subtree pull --prefix "${prefixes[$remote]}" --squash "${urls[$remote]}" master
|
|
||||||
done
|
|
||||||
|
|
||||||
for remote in "${!urls[@]}"; do
|
|
||||||
git_subtree "push" "${prefixes[$remote]}" "$remote" "${urls[$remote]}"
|
|
||||||
done
|
|
||||||
|
|
||||||
git push origin master
|
|
||||||
git status
|
|
||||||
|
|
||||||
trigger-downstream:
|
|
||||||
stage: trigger
|
|
||||||
##
|
|
||||||
## Only run if on the master branch and the variable GROUP is set
|
|
||||||
##
|
|
||||||
## change this to
|
|
||||||
## only:
|
|
||||||
## - master@$GROUP/iphreeqc
|
|
||||||
## and set GROUP to coupled before merge
|
|
||||||
only:
|
|
||||||
refs:
|
|
||||||
- master
|
|
||||||
variables:
|
|
||||||
- $GROUP
|
|
||||||
|
|
||||||
## Downstream Projects
|
|
||||||
## triggers and ids are stored at the group level
|
|
||||||
## webmod https://code.chs.usgs.gov/coupled/webmod
|
|
||||||
script:
|
|
||||||
- echo triggering webmod
|
|
||||||
- curl -X POST -F token=${WEBMOD_TRIGGER} -F ref=master https://code.chs.usgs.gov/api/v4/projects/${WEBMOD_ID}/trigger/pipeline
|
|
||||||
- sleep 180
|
|
||||||
|
|
||||||
## Upstream Projects
|
|
||||||
## iphreeqc-src https://code.chs.usgs.gov/coupled/subtrees/iphreeqc-src
|
|
||||||
## phreeqc-commanuscript-cgfinal-examples-c https://code.chs.usgs.gov/coupled/subtrees/phreeqc-commanuscript-cgfinal-examples-c
|
|
||||||
## phreeqc-commanuscript-cgfinal-examples-com https://code.chs.usgs.gov/coupled/subtrees/phreeqc-commanuscript-cgfinal-examples-com
|
|
||||||
## phreeqc-COMManuscript-CGfinal-examples-fortran https://code.chs.usgs.gov/coupled/subtrees/phreeqc-COMManuscript-CGfinal-examples-fortran
|
|
||||||
## phreeqc3-database https://code.chs.usgs.gov/coupled/subtrees/phreeqc3-database
|
|
||||||
## phreeqc3-doc https://code.chs.usgs.gov/coupled/subtrees/phreeqc3-doc
|
|
||||||
## phreeqc3-examples https://code.chs.usgs.gov/coupled/subtrees/phreeqc3-examples
|
|
||||||
@ -220,6 +220,7 @@ target_sources(IPhreeqc
|
|||||||
src/phreeqcpp/UserPunch.cpp
|
src/phreeqcpp/UserPunch.cpp
|
||||||
src/phreeqcpp/UserPunch.h
|
src/phreeqcpp/UserPunch.h
|
||||||
src/phreeqcpp/utilities.cpp
|
src/phreeqcpp/utilities.cpp
|
||||||
|
src/phreeqcpp/GFZ.cpp
|
||||||
src/thread.h
|
src/thread.h
|
||||||
src/Var.c
|
src/Var.c
|
||||||
src/Var.h
|
src/Var.h
|
||||||
@ -392,9 +393,6 @@ if (BUILD_CLR_LIBS)
|
|||||||
|
|
||||||
endif()
|
endif()
|
||||||
|
|
||||||
add_subdirectory(poet)
|
|
||||||
add_subdirectory(app)
|
|
||||||
|
|
||||||
include (CTest)
|
include (CTest)
|
||||||
|
|
||||||
if (STANDALONE_BUILD)
|
if (STANDALONE_BUILD)
|
||||||
@ -459,6 +457,9 @@ if (STANDALONE_BUILD)
|
|||||||
endif()
|
endif()
|
||||||
endif()
|
endif()
|
||||||
|
|
||||||
|
add_subdirectory(poet)
|
||||||
|
# add_subdirectory(app)
|
||||||
|
|
||||||
# get_cmake_property(_variableNames VARIABLES)
|
# get_cmake_property(_variableNames VARIABLES)
|
||||||
# list (SORT _variableNames)
|
# list (SORT _variableNames)
|
||||||
# foreach (_variableName ${_variableNames})
|
# foreach (_variableName ${_variableNames})
|
||||||
|
|||||||
@ -5,8 +5,11 @@ add_library(IPhreeqcPOET ${PQC_POET_SRC})
|
|||||||
target_link_libraries(IPhreeqcPOET PUBLIC IPhreeqc)
|
target_link_libraries(IPhreeqcPOET PUBLIC IPhreeqc)
|
||||||
target_include_directories(IPhreeqcPOET PUBLIC
|
target_include_directories(IPhreeqcPOET PUBLIC
|
||||||
$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
|
$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
|
||||||
$<INSTALL_INTERFACE:include>
|
$<INSTALL_INTERFACE:include>)
|
||||||
PRIVATE src)
|
|
||||||
|
|
||||||
# Set C++23 standard
|
# Set C++20 standard
|
||||||
target_compile_features(IPhreeqcPOET PUBLIC cxx_std_20)
|
target_compile_features(IPhreeqcPOET PUBLIC cxx_std_20)
|
||||||
|
|
||||||
|
if (BUILD_TESTING AND STANDALONE_BUILD)
|
||||||
|
add_subdirectory(test)
|
||||||
|
endif()
|
||||||
@ -1,20 +0,0 @@
|
|||||||
#pragma once
|
|
||||||
|
|
||||||
#include <string>
|
|
||||||
#include <vector>
|
|
||||||
|
|
||||||
struct POETInitCell {
|
|
||||||
std::vector<std::string> solutions;
|
|
||||||
std::vector<std::string> solution_primaries;
|
|
||||||
std::vector<std::string> exchanger;
|
|
||||||
std::vector<std::string> kinetics;
|
|
||||||
std::vector<std::string> equilibrium;
|
|
||||||
std::vector<std::string> surface_comps;
|
|
||||||
std::vector<std::string> surface_charges;
|
|
||||||
};
|
|
||||||
|
|
||||||
struct POETConfig {
|
|
||||||
std::string database;
|
|
||||||
std::string input_script;
|
|
||||||
POETInitCell cell;
|
|
||||||
};
|
|
||||||
@ -1,16 +1,15 @@
|
|||||||
#pragma once
|
#pragma once
|
||||||
|
|
||||||
#include "POETInit.hpp"
|
#include "PhreeqcMatrix.hpp"
|
||||||
|
|
||||||
#include <memory>
|
#include <memory>
|
||||||
#include <vector>
|
#include <vector>
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* @brief Class for running Phreeqc wrappped in POET
|
* @brief Class for running Phreeqc wrappped in POET
|
||||||
*
|
*
|
||||||
* Direct interface to Phreeqc, without utilizing *_MODIFY keywords/scripts
|
* Direct interface to Phreeqc, without utilizing Phreeqc's internal parser
|
||||||
to
|
* to set/get values. To start a simulation, the interpreter is used internally.
|
||||||
* set new values. Use with already initialized Phreeqc config.
|
* No need for the user to directly interact with the interpreter.
|
||||||
*
|
*
|
||||||
*/
|
*/
|
||||||
class PhreeqcEngine {
|
class PhreeqcEngine {
|
||||||
@ -19,11 +18,13 @@ public:
|
|||||||
* @brief Construct a new Phreeqc Engine object
|
* @brief Construct a new Phreeqc Engine object
|
||||||
*
|
*
|
||||||
* Construct a new Phreeqc Engine object by previously initialized
|
* Construct a new Phreeqc Engine object by previously initialized
|
||||||
POETConfig.
|
* PhreeqcMatrix.
|
||||||
*
|
*
|
||||||
* @param config Holds the configuration for the Phreeqc engine.
|
* @param pqc_mat PhreeqcMatrix initialized with the *full* Phreeqc script and
|
||||||
|
* database.
|
||||||
|
* @param cell_id ID of the cell (user id from Phreeqc script) to simulate
|
||||||
*/
|
*/
|
||||||
PhreeqcEngine(const POETConfig &config);
|
PhreeqcEngine(const PhreeqcMatrix &pqc_mat, const int cell_id);
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* @brief Destroy the Phreeqc Engine object
|
* @brief Destroy the Phreeqc Engine object
|
||||||
@ -35,8 +36,8 @@ public:
|
|||||||
* @brief Siimulate a cell for a given time step
|
* @brief Siimulate a cell for a given time step
|
||||||
*
|
*
|
||||||
* @param cell_values Vector containing the input values for the cell
|
* @param cell_values Vector containing the input values for the cell
|
||||||
* (**including the ID**). Output values are written back in place to this
|
* (*including the ID*). Output values are written back in place to this
|
||||||
* vector!
|
* vector.
|
||||||
* @param time_step Time step to simulate in seconds
|
* @param time_step Time step to simulate in seconds
|
||||||
*/
|
*/
|
||||||
void runCell(std::vector<double> &cell_values, double time_step);
|
void runCell(std::vector<double> &cell_values, double time_step);
|
||||||
|
|||||||
@ -1,247 +0,0 @@
|
|||||||
#pragma once
|
|
||||||
|
|
||||||
#include <IPhreeqc.hpp>
|
|
||||||
|
|
||||||
#include <Phreeqc.h>
|
|
||||||
#include <array>
|
|
||||||
#include <cstddef>
|
|
||||||
#include <cstdint>
|
|
||||||
#include <map>
|
|
||||||
#include <string>
|
|
||||||
#include <utility>
|
|
||||||
#include <vector>
|
|
||||||
|
|
||||||
#include "Exchange.h"
|
|
||||||
#include "PPassemblage.h"
|
|
||||||
#include "Solution.h"
|
|
||||||
#include "Surface.h"
|
|
||||||
#include "cxxKinetics.h"
|
|
||||||
|
|
||||||
enum { POET_SOL = 0, POET_EXCH, POET_KIN, POET_EQUIL, POET_SURF };
|
|
||||||
|
|
||||||
/**
|
|
||||||
* @brief Class for initializing Phreeqc and create input for POET
|
|
||||||
*
|
|
||||||
* This class is used to initialize Phreeqc by reading a database and an input
|
|
||||||
* Phreeqc script. It then extracts the 'plain' initial values for each cell and
|
|
||||||
* module, which are then used for further processing in POET.
|
|
||||||
*/
|
|
||||||
class PhreeqcInit : public IPhreeqc {
|
|
||||||
public:
|
|
||||||
/**
|
|
||||||
* @brief Construct a new Phreeqc Init object
|
|
||||||
*
|
|
||||||
* @param database String containing the content of the database file
|
|
||||||
* @param input_script String containing the content of the input script
|
|
||||||
*/
|
|
||||||
PhreeqcInit(const std::string &database, const std::string &input_script);
|
|
||||||
|
|
||||||
/**
|
|
||||||
* Hold column names, indexes and values for each cell defined in the input
|
|
||||||
* script
|
|
||||||
*/
|
|
||||||
struct PhreeqcMat {
|
|
||||||
std::vector<std::string> names;
|
|
||||||
std::vector<int> ids;
|
|
||||||
std::vector<std::vector<double>> values;
|
|
||||||
};
|
|
||||||
|
|
||||||
/**
|
|
||||||
* @brief Get the PhreeqcMatrix
|
|
||||||
*
|
|
||||||
* @return PhreeqcMat Value matrix of the parsed Phreeqc input script
|
|
||||||
*/
|
|
||||||
PhreeqcMat getPhreeqcMat() const { return pqc_mat; };
|
|
||||||
|
|
||||||
/**
|
|
||||||
* @brief Get the *_RAW input scripts for each cell
|
|
||||||
*
|
|
||||||
* @return std::map<int, std::string> Contains the raw input scripts for each
|
|
||||||
* cell
|
|
||||||
*/
|
|
||||||
std::map<int, std::string> raw_dumps();
|
|
||||||
|
|
||||||
using essential_names = std::array<std::vector<std::string>, 5>;
|
|
||||||
using ModulesArray = std::array<std::uint32_t, 5>;
|
|
||||||
|
|
||||||
/**
|
|
||||||
* @brief Get the actual reactant sizes for given cells
|
|
||||||
*
|
|
||||||
* @param cell_ids Vector of cell ids to return reactant sizes for
|
|
||||||
* @return ModulesArray Combined sizes of all cells for each module
|
|
||||||
*/
|
|
||||||
ModulesArray getModuleSizes(const std::vector<int> &cell_ids);
|
|
||||||
|
|
||||||
/**
|
|
||||||
* @brief Get solution primaries found by PhreeqcInit
|
|
||||||
*
|
|
||||||
* @return std::vector<std::string> Vector containing only the solution
|
|
||||||
* primary names of the solution defined by the input script
|
|
||||||
*/
|
|
||||||
std::vector<std::string> getSolutionPrimaries() {
|
|
||||||
return std::vector<std::string>(this->surface_primaries.begin(),
|
|
||||||
this->surface_primaries.end());
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
|
||||||
* @brief Get the solution names for a given cell
|
|
||||||
*
|
|
||||||
* @param cell_id ID of the cell to get the solution names for
|
|
||||||
* @return std::vector<std::string> Whole vector of solution names for the
|
|
||||||
* cell
|
|
||||||
*/
|
|
||||||
std::vector<std::string> getSolutionNames(int cell_id) {
|
|
||||||
return this->raw_initials[cell_id].first[POET_SOL];
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
|
||||||
* @brief Get the exchange names for a given cell
|
|
||||||
*
|
|
||||||
* @param cell_id ID of the cell to get the exchange names for
|
|
||||||
* @return std::vector<std::string> Whole vector of exchange names for the
|
|
||||||
* cell. Empty if no exchange is defined.
|
|
||||||
*/
|
|
||||||
std::vector<std::string> getExchanger(int id) {
|
|
||||||
if (contains_id(this->exchanger, id)) {
|
|
||||||
return this->exchanger[id];
|
|
||||||
}
|
|
||||||
|
|
||||||
return {};
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
|
||||||
* @brief Get the kinetics names for a given cell
|
|
||||||
*
|
|
||||||
* @param cell_id ID of the cell to get the kinetics names for
|
|
||||||
* @return std::vector<std::string> Whole vector of kinetics names for the
|
|
||||||
* cell. Empty if no kinetics are defined.
|
|
||||||
*/
|
|
||||||
std::vector<std::string> getKineticsNames(int id) {
|
|
||||||
if (contains_id(this->kinetics, id)) {
|
|
||||||
return this->kinetics[id];
|
|
||||||
}
|
|
||||||
|
|
||||||
return {};
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
|
||||||
* @brief Get the equilibrium names for a given cell
|
|
||||||
*
|
|
||||||
* @param cell_id ID of the cell to get the equilibrium names for
|
|
||||||
* @return std::vector<std::string> Whole vector of equilibrium names for the
|
|
||||||
* cell. Empty if no equilibrium is defined.
|
|
||||||
*/
|
|
||||||
std::vector<std::string> getEquilibriumNames(int id) {
|
|
||||||
if (contains_id(this->equilibrium, id)) {
|
|
||||||
return this->equilibrium[id];
|
|
||||||
}
|
|
||||||
|
|
||||||
return {};
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
|
||||||
* @brief Get the surface component names for a given cell
|
|
||||||
*
|
|
||||||
* @param cell_id ID of the cell to get the surface component names for
|
|
||||||
* @return std::vector<std::string> Whole vector of surface component names
|
|
||||||
* for the cell. Empty if no surface is defined.
|
|
||||||
*/
|
|
||||||
std::vector<std::string> getSurfaceCompNames(int id) {
|
|
||||||
if (contains_id(this->surface_comps, id)) {
|
|
||||||
return this->surface_comps[id];
|
|
||||||
}
|
|
||||||
|
|
||||||
return {};
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
|
||||||
* @brief Get the surface charge names for a given cell
|
|
||||||
*
|
|
||||||
* @param cell_id ID of the cell to get the surface charge names for
|
|
||||||
* @return std::vector<std::string> Whole vector of surface charge names for
|
|
||||||
* the cell. Empty if no surface is defined.
|
|
||||||
*/
|
|
||||||
std::vector<std::string> getSurfaceChargeNames(int id) {
|
|
||||||
if (contains_id(this->surface_charge, id)) {
|
|
||||||
return this->surface_charge[id];
|
|
||||||
}
|
|
||||||
|
|
||||||
return {};
|
|
||||||
}
|
|
||||||
|
|
||||||
private:
|
|
||||||
// required only for simulation
|
|
||||||
essential_names initial_names;
|
|
||||||
|
|
||||||
void parseInitValues();
|
|
||||||
void parseInit();
|
|
||||||
|
|
||||||
std::vector<std::string> dump_solution_names(int cell_number);
|
|
||||||
void dump_reactant_names(int cell_number,
|
|
||||||
const std::vector<std::string> union_sol_names,
|
|
||||||
essential_names &names);
|
|
||||||
|
|
||||||
std::vector<std::string>
|
|
||||||
find_all_valence_states(const std::vector<std::string> &solution_names,
|
|
||||||
const std::size_t offset);
|
|
||||||
|
|
||||||
cxxSolution *Get_solution(std::size_t n) {
|
|
||||||
return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_solution_map(), n);
|
|
||||||
}
|
|
||||||
|
|
||||||
cxxExchange *Get_exchange(std::size_t n) {
|
|
||||||
return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_exchange_map(), n);
|
|
||||||
}
|
|
||||||
|
|
||||||
cxxKinetics *Get_kinetic(std::size_t n) {
|
|
||||||
return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_kinetics_map(), n);
|
|
||||||
}
|
|
||||||
|
|
||||||
cxxPPassemblage *Get_equilibrium(std::size_t n) {
|
|
||||||
return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_pp_assemblage_map(),
|
|
||||||
n);
|
|
||||||
}
|
|
||||||
|
|
||||||
cxxSurface *Get_surface(std::size_t n) {
|
|
||||||
return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_surface_map(), n);
|
|
||||||
}
|
|
||||||
|
|
||||||
PhreeqcMat pqc_mat;
|
|
||||||
|
|
||||||
PhreeqcMat buildPhreeqcMat();
|
|
||||||
|
|
||||||
std::map<int, std::vector<std::string>> exchanger;
|
|
||||||
std::map<int, std::vector<std::string>> kinetics;
|
|
||||||
std::map<int, std::vector<std::string>> equilibrium;
|
|
||||||
std::map<int, std::vector<std::string>> surface_comps;
|
|
||||||
std::map<int, std::vector<std::string>> surface_charge;
|
|
||||||
|
|
||||||
bool contains_id(const std::map<int, std::vector<std::string>> &map,
|
|
||||||
int id) const {
|
|
||||||
#if __cplusplus >= 202002L
|
|
||||||
return map.contains(id);
|
|
||||||
#else
|
|
||||||
return map.find(id) != map.end();
|
|
||||||
#endif
|
|
||||||
}
|
|
||||||
|
|
||||||
std::set<std::string> surface_primaries;
|
|
||||||
|
|
||||||
using RawMap = std::map<int, std::pair<essential_names, std::vector<double>>>;
|
|
||||||
|
|
||||||
essential_names union_raws(const RawMap &raws);
|
|
||||||
|
|
||||||
std::vector<std::vector<double>>
|
|
||||||
conc_from_essentials(const RawMap &raws, const essential_names &names);
|
|
||||||
|
|
||||||
// void valuesFromModule(const std::string &module_name, int cell_number,
|
|
||||||
// essential_names &names, std::vector<double> &values);
|
|
||||||
|
|
||||||
std::string subExchangeName(std::string name);
|
|
||||||
|
|
||||||
std::vector<double>
|
|
||||||
get_essential_values_init(std::size_t cell_number,
|
|
||||||
const std::vector<std::string> &order);
|
|
||||||
|
|
||||||
RawMap raw_initials;
|
|
||||||
};
|
|
||||||
274
poet/include/PhreeqcMatrix.hpp
Normal file
274
poet/include/PhreeqcMatrix.hpp
Normal file
@ -0,0 +1,274 @@
|
|||||||
|
#pragma once
|
||||||
|
|
||||||
|
#include <array>
|
||||||
|
#include <cstddef>
|
||||||
|
#include <map>
|
||||||
|
#include <memory>
|
||||||
|
#include <set>
|
||||||
|
#include <string>
|
||||||
|
#include <vector>
|
||||||
|
|
||||||
|
#include "IPhreeqc.hpp"
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @brief Class for storing information from Phreeqc
|
||||||
|
*
|
||||||
|
* PhreeqcMatrix is used as a container for storing **essential** information
|
||||||
|
* from Phreeqc in a C++ data structure. The usage of Phreeqc's interpreter is
|
||||||
|
* minimized. Thus, values are written directly from/to Phreeqc's internal data
|
||||||
|
* structures, minimizing the overhead of parsing and eliminate floating point
|
||||||
|
* errors due to conversion.
|
||||||
|
*
|
||||||
|
* The class is also used to initialize the PhreeqcEngine class.
|
||||||
|
*/
|
||||||
|
class PhreeqcMatrix {
|
||||||
|
public:
|
||||||
|
/**
|
||||||
|
* @brief Construct a new Phreeqc Matrix object
|
||||||
|
*
|
||||||
|
* Default constructor. Does nothing. Used only for assignment operator.
|
||||||
|
*/
|
||||||
|
PhreeqcMatrix() = default;
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @brief Construct a new Phreeqc Matrix object
|
||||||
|
*
|
||||||
|
* Construct a new Phreeqc Matrix object by reading the database and input
|
||||||
|
* script already present as a string.
|
||||||
|
*
|
||||||
|
* @param database
|
||||||
|
* @param input_script
|
||||||
|
*/
|
||||||
|
PhreeqcMatrix(const std::string &database, const std::string &input_script);
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @brief Construct a new Phreeqc Matrix object
|
||||||
|
*
|
||||||
|
* Copy constructor. The interal used Phreeqc instance is reused by the new
|
||||||
|
* object!
|
||||||
|
* @param other
|
||||||
|
*/
|
||||||
|
PhreeqcMatrix(const PhreeqcMatrix &other);
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @brief Construct a new Phreeqc Matrix object
|
||||||
|
*
|
||||||
|
* Copy constructor. The interal used Phreeqc instance is reused by the new
|
||||||
|
* object!
|
||||||
|
* @param other
|
||||||
|
*/
|
||||||
|
PhreeqcMatrix(PhreeqcMatrix &&other);
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @brief Assignment operator
|
||||||
|
*
|
||||||
|
* The interal used Phreeqc instance is reused by the new object!
|
||||||
|
* @param other
|
||||||
|
* @return PhreeqcMatrix&
|
||||||
|
*/
|
||||||
|
PhreeqcMatrix &operator=(const PhreeqcMatrix &other);
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @brief Assignment operator
|
||||||
|
*
|
||||||
|
* The interal used Phreeqc instance is reused by the new object!
|
||||||
|
* @param other
|
||||||
|
* @return PhreeqcMatrix&
|
||||||
|
*/
|
||||||
|
PhreeqcMatrix &operator=(PhreeqcMatrix &&other);
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @brief Access the value of a given cell by name.
|
||||||
|
*
|
||||||
|
* @param cell_id ID of the cell (user id from Phreeqc script) to get the
|
||||||
|
* value from.
|
||||||
|
* @param name Name of the component to get the value from.
|
||||||
|
* @return double The stored value.
|
||||||
|
* @throw std::runtime_error if the component is not found.
|
||||||
|
*/
|
||||||
|
double operator()(int cell_id, const std::string &name) const;
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @brief Subset the PhreeqcMatrix to given cell IDs.
|
||||||
|
*
|
||||||
|
* With a given set of cell IDs, a new PhreeqcMatrix is created containing
|
||||||
|
* only the given cell IDs. All entries, which values refer to NaN, are
|
||||||
|
* removed.
|
||||||
|
*
|
||||||
|
* @param indices Cell IDs to subset the PhreeqcMatrix to.
|
||||||
|
* @return PhreeqcMatrix A new PhreeqcMatrix containing only the given cell
|
||||||
|
* IDs.
|
||||||
|
*/
|
||||||
|
PhreeqcMatrix subset(const std::vector<int> &indices) const;
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @brief Erase the given cell IDs from the PhreeqcMatrix.
|
||||||
|
*
|
||||||
|
* With a given set of cell IDs, the PhreeqcMatrix is modified to contain only
|
||||||
|
* the cell IDs not in the given set. All entries, which values refer to NaN,
|
||||||
|
* are removed.
|
||||||
|
*
|
||||||
|
* @param indices Cell IDs to erase from the PhreeqcMatrix.
|
||||||
|
* @return PhreeqcMatrix A new PhreeqcMatrix containing only the cell IDs not
|
||||||
|
* in the given set.
|
||||||
|
*/
|
||||||
|
PhreeqcMatrix erase(const std::vector<int> &indices) const;
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @brief Type of vector export
|
||||||
|
*
|
||||||
|
*/
|
||||||
|
enum class VectorExportType { COLUMN_MAJOR, ROW_MAJOR };
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @brief Struct holding a format of the exported data
|
||||||
|
*
|
||||||
|
*/
|
||||||
|
struct STLExport {
|
||||||
|
std::vector<std::string> names;
|
||||||
|
std::vector<double> values;
|
||||||
|
};
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @brief Export the internal data to consecutive vectors.
|
||||||
|
*
|
||||||
|
* @param type Type of the order of the exported data
|
||||||
|
* @param include_id Whether to include a column with the cell IDs or not
|
||||||
|
* @return STLExport Exported data
|
||||||
|
*/
|
||||||
|
STLExport get(VectorExportType type = VectorExportType::ROW_MAJOR,
|
||||||
|
bool include_id = true) const;
|
||||||
|
|
||||||
|
enum class PhreeqcComponent {
|
||||||
|
SOLUTION = 0,
|
||||||
|
EXCHANGE,
|
||||||
|
KINETIC,
|
||||||
|
EQUILIBRIUM,
|
||||||
|
SURFACE_COMPS
|
||||||
|
};
|
||||||
|
|
||||||
|
struct element {
|
||||||
|
std::string name;
|
||||||
|
PhreeqcComponent type;
|
||||||
|
double value;
|
||||||
|
};
|
||||||
|
|
||||||
|
struct base_names {
|
||||||
|
enum class Components {
|
||||||
|
EXCHANGER = static_cast<int>(PhreeqcComponent::EXCHANGE),
|
||||||
|
KINETICS,
|
||||||
|
EQUILIBRIUM,
|
||||||
|
SURACE_COMP,
|
||||||
|
SURFACE_CHARGE
|
||||||
|
} type;
|
||||||
|
|
||||||
|
std::string name;
|
||||||
|
};
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @brief Get all found solution names
|
||||||
|
*
|
||||||
|
* @return std::vector<std::string> Vector containing all solution names.
|
||||||
|
*/
|
||||||
|
std::vector<std::string>
|
||||||
|
getSolutionNames(bool include_h_o_charge = false) const;
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @brief Get solution total names of all found solutions (excluding H, O,
|
||||||
|
* Charge)
|
||||||
|
*
|
||||||
|
* @return std::vector<std::string> Names of all found solutions (excluding H,
|
||||||
|
* O, Charge)
|
||||||
|
*/
|
||||||
|
std::vector<std::string> getSolutionPrimaries() const;
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @brief Get the exchange names for a given cell
|
||||||
|
*
|
||||||
|
* @param cell_id ID of the cell to get the exchange names for
|
||||||
|
* @return std::vector<std::string> Whole vector of exchange names for the
|
||||||
|
* cell. Empty if no exchange is defined.
|
||||||
|
*/
|
||||||
|
std::vector<std::string> getExchanger(int cell_id) const;
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @brief Get the kinetics names for a given cell
|
||||||
|
*
|
||||||
|
* @param cell_id ID of the cell to get the kinetics names for
|
||||||
|
* @return std::vector<std::string> Whole vector of kinetics names for the
|
||||||
|
* cell. Empty if no kinetics are defined.
|
||||||
|
*/
|
||||||
|
std::vector<std::string> getKineticsNames(int cell_id) const;
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @brief Get the equilibrium names for a given cell
|
||||||
|
*
|
||||||
|
* @param cell_id ID of the cell to get the equilibrium names for
|
||||||
|
* @return std::vector<std::string> Whole vector of equilibrium names for the
|
||||||
|
* cell. Empty if no equilibrium is defined.
|
||||||
|
*/
|
||||||
|
std::vector<std::string> getEquilibriumNames(int cell_id) const;
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @brief Get the surface component names for a given cell
|
||||||
|
*
|
||||||
|
* @param cell_id ID of the cell to get the surface component names for
|
||||||
|
* @return std::vector<std::string> Whole vector of surface component names
|
||||||
|
* for the cell. Empty if no surface is defined.
|
||||||
|
*/
|
||||||
|
std::vector<std::string> getSurfaceCompNames(int cell_id) const;
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @brief Get the surface charge names for a given cell
|
||||||
|
*
|
||||||
|
* @param cell_id ID of the cell to get the surface charge names for
|
||||||
|
* @return std::vector<std::string> Whole vector of surface charge names for
|
||||||
|
* the cell. Empty if no surface is defined.
|
||||||
|
*/
|
||||||
|
std::vector<std::string> getSurfaceChargeNames(int cell_id) const;
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @brief Get all cell IDs stored in the PhreeqcMatrix.
|
||||||
|
*
|
||||||
|
* @return std::vector<int> IDs of all cells stored in the PhreeqcMatrix.
|
||||||
|
*/
|
||||||
|
std::vector<int> getIds() const;
|
||||||
|
|
||||||
|
// std::array<std::size_t, 5> getComponentCount(int cell_id) const;
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @brief Dump all cells into a **DUMP** format of Phreeqc.
|
||||||
|
*
|
||||||
|
* @return std::map<int, std::string> Map containing the cell ID as key and
|
||||||
|
* the exported DUMP string as value.
|
||||||
|
*/
|
||||||
|
std::map<int, std::string> getDumpStringsPQI() const;
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @brief Get the **DUMP** string for a given cell.
|
||||||
|
*
|
||||||
|
* @param cell_id Cell ID to get the **DUMP** string for.
|
||||||
|
* @return std::string Phreeqc **DUMP** string for the given cell.
|
||||||
|
*/
|
||||||
|
std::string getDumpStringsPQI(int cell_id) const;
|
||||||
|
|
||||||
|
/**
|
||||||
|
* @brief Get the Database used to initialize the PhreeqcMatrix.
|
||||||
|
*
|
||||||
|
* @return std::string Database string.
|
||||||
|
*/
|
||||||
|
std::string getDatabase() const;
|
||||||
|
|
||||||
|
private:
|
||||||
|
std::map<int, std::vector<element>> _m_map;
|
||||||
|
std::map<int, std::vector<base_names>> _m_internal_names;
|
||||||
|
|
||||||
|
std::set<std::string> _m_surface_primaries;
|
||||||
|
|
||||||
|
void initialize();
|
||||||
|
|
||||||
|
void remove_NaNs();
|
||||||
|
|
||||||
|
std::shared_ptr<IPhreeqc> _m_pqc;
|
||||||
|
std::string _m_database;
|
||||||
|
};
|
||||||
@ -1,12 +1,16 @@
|
|||||||
#include "PhreeqcEngine.hpp"
|
#include "PhreeqcEngine.hpp"
|
||||||
#include <cstddef>
|
#include <cstddef>
|
||||||
#include <iomanip>
|
#include <iomanip>
|
||||||
|
#include <regex>
|
||||||
#include <span>
|
#include <span>
|
||||||
#include <sstream>
|
#include <sstream>
|
||||||
|
|
||||||
#include <IPhreeqc.hpp>
|
#include <IPhreeqc.hpp>
|
||||||
#include <Phreeqc.h>
|
#include <Phreeqc.h>
|
||||||
|
#include <string>
|
||||||
|
#include <vector>
|
||||||
|
|
||||||
|
#include "PhreeqcMatrix.hpp"
|
||||||
#include "Wrapper/EquilibriumWrapper.hpp"
|
#include "Wrapper/EquilibriumWrapper.hpp"
|
||||||
#include "Wrapper/ExchangeWrapper.hpp"
|
#include "Wrapper/ExchangeWrapper.hpp"
|
||||||
#include "Wrapper/KineticWrapper.hpp"
|
#include "Wrapper/KineticWrapper.hpp"
|
||||||
@ -14,8 +18,7 @@
|
|||||||
#include "Wrapper/SurfaceWrapper.hpp"
|
#include "Wrapper/SurfaceWrapper.hpp"
|
||||||
class PhreeqcEngine::Impl : public IPhreeqc {
|
class PhreeqcEngine::Impl : public IPhreeqc {
|
||||||
public:
|
public:
|
||||||
Impl(const POETConfig &config);
|
Impl(const PhreeqcMatrix &pqc_mat, const int cell_id);
|
||||||
void init_wrappers(const POETInitCell &cell);
|
|
||||||
void run(double time_step);
|
void run(double time_step);
|
||||||
|
|
||||||
cxxSolution *Get_solution(std::size_t n) {
|
cxxSolution *Get_solution(std::size_t n) {
|
||||||
@ -53,20 +56,48 @@ public:
|
|||||||
bool has_kinetics = false;
|
bool has_kinetics = false;
|
||||||
bool has_equilibrium = false;
|
bool has_equilibrium = false;
|
||||||
bool has_surface = false;
|
bool has_surface = false;
|
||||||
|
struct InitCell {
|
||||||
|
std::vector<std::string> solutions;
|
||||||
|
std::vector<std::string> exchanger;
|
||||||
|
std::vector<std::string> kinetics;
|
||||||
|
std::vector<std::string> equilibrium;
|
||||||
|
std::vector<std::string> surface_comps;
|
||||||
|
std::vector<std::string> surface_charges;
|
||||||
|
std::vector<std::string> solution_primaries;
|
||||||
|
};
|
||||||
|
void init_wrappers(const InitCell &cell);
|
||||||
};
|
};
|
||||||
|
|
||||||
PhreeqcEngine::PhreeqcEngine(const POETConfig &config)
|
PhreeqcEngine::PhreeqcEngine(const PhreeqcMatrix &pqc_mat, const int cell_id)
|
||||||
: impl(std::make_unique<Impl>(config)) {}
|
: impl(std::make_unique<Impl>(pqc_mat, cell_id)) {}
|
||||||
|
|
||||||
PhreeqcEngine::Impl::Impl(const POETConfig &config) {
|
|
||||||
this->LoadDatabaseString(config.database.c_str());
|
|
||||||
this->RunString(config.input_script.c_str());
|
|
||||||
|
|
||||||
this->init_wrappers(config.cell);
|
|
||||||
}
|
|
||||||
|
|
||||||
PhreeqcEngine::~PhreeqcEngine() = default;
|
PhreeqcEngine::~PhreeqcEngine() = default;
|
||||||
|
|
||||||
|
static inline std::string
|
||||||
|
replaceRawKeywordID(const std::string &raw_dump_string) {
|
||||||
|
std::regex re(R"((RAW\s+)(\d+))");
|
||||||
|
return std::regex_replace(raw_dump_string, re, "RAW 1");
|
||||||
|
}
|
||||||
|
|
||||||
|
PhreeqcEngine::Impl::Impl(const PhreeqcMatrix &pqc_mat, const int cell_id) {
|
||||||
|
this->LoadDatabaseString(pqc_mat.getDatabase().c_str());
|
||||||
|
|
||||||
|
const std::string pqc_string =
|
||||||
|
replaceRawKeywordID(pqc_mat.getDumpStringsPQI(cell_id));
|
||||||
|
|
||||||
|
this->RunString(pqc_string.c_str());
|
||||||
|
|
||||||
|
InitCell cell = {pqc_mat.getSolutionNames(),
|
||||||
|
pqc_mat.getExchanger(cell_id),
|
||||||
|
pqc_mat.getKineticsNames(cell_id),
|
||||||
|
pqc_mat.getEquilibriumNames(cell_id),
|
||||||
|
pqc_mat.getSurfaceCompNames(cell_id),
|
||||||
|
pqc_mat.getSurfaceChargeNames(cell_id),
|
||||||
|
pqc_mat.getSolutionPrimaries()};
|
||||||
|
|
||||||
|
this->init_wrappers(cell);
|
||||||
|
}
|
||||||
|
|
||||||
void PhreeqcEngine::runCell(std::vector<double> &cell_values,
|
void PhreeqcEngine::runCell(std::vector<double> &cell_values,
|
||||||
double time_step) {
|
double time_step) {
|
||||||
// skip ID
|
// skip ID
|
||||||
@ -86,7 +117,7 @@ void PhreeqcEngine::Impl::run(double time_step) {
|
|||||||
this->RunString(runs_string.c_str());
|
this->RunString(runs_string.c_str());
|
||||||
}
|
}
|
||||||
|
|
||||||
void PhreeqcEngine::Impl::init_wrappers(const POETInitCell &cell) {
|
void PhreeqcEngine::Impl::init_wrappers(const InitCell &cell) {
|
||||||
|
|
||||||
// Solutions
|
// Solutions
|
||||||
this->solutionWrapperPtr =
|
this->solutionWrapperPtr =
|
||||||
|
|||||||
@ -1,204 +0,0 @@
|
|||||||
#include <PhreeqcInit.hpp>
|
|
||||||
|
|
||||||
#include "Solution.h"
|
|
||||||
|
|
||||||
#include <algorithm>
|
|
||||||
#include <cstddef>
|
|
||||||
#include <tuple>
|
|
||||||
#include <vector>
|
|
||||||
|
|
||||||
static inline std::vector<std::string>
|
|
||||||
unionStringVectors(const std::vector<std::string> &a,
|
|
||||||
const std::vector<std::string> &b) {
|
|
||||||
std::vector<std::string> result;
|
|
||||||
|
|
||||||
std::set_union(a.begin(), a.end(), b.begin(), b.end(),
|
|
||||||
std::back_inserter(result));
|
|
||||||
|
|
||||||
return result;
|
|
||||||
}
|
|
||||||
|
|
||||||
static std::vector<double>
|
|
||||||
createConcVector(const std::vector<std::string> &conc_names,
|
|
||||||
const std::vector<double>::const_iterator &conc_it,
|
|
||||||
const std::vector<std::string> &new_names) {
|
|
||||||
std::vector<double> conc_vec;
|
|
||||||
|
|
||||||
for (const auto &i : new_names) {
|
|
||||||
auto it = std::find(conc_names.begin(), conc_names.end(), i);
|
|
||||||
|
|
||||||
auto conc_index = std::distance(conc_names.begin(), it);
|
|
||||||
|
|
||||||
if (it != conc_names.end()) {
|
|
||||||
conc_vec.push_back(*(conc_it + conc_index));
|
|
||||||
} else {
|
|
||||||
conc_vec.push_back(std::numeric_limits<double>::quiet_NaN());
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
return conc_vec;
|
|
||||||
}
|
|
||||||
|
|
||||||
PhreeqcInit::PhreeqcInit(const std::string &database,
|
|
||||||
const std::string &input_script)
|
|
||||||
: IPhreeqc() {
|
|
||||||
this->LoadDatabaseString(database.c_str());
|
|
||||||
this->RunString(input_script.c_str());
|
|
||||||
|
|
||||||
this->parseInit();
|
|
||||||
pqc_mat = this->buildPhreeqcMat();
|
|
||||||
}
|
|
||||||
|
|
||||||
void PhreeqcInit::parseInit() {
|
|
||||||
|
|
||||||
std::map<int, essential_names> init_values;
|
|
||||||
|
|
||||||
for (const auto &[id, val] : this->PhreeqcPtr->Get_Rxn_solution_map()) {
|
|
||||||
// A key less than zero indicates an internal solution
|
|
||||||
if (id < 0) {
|
|
||||||
continue;
|
|
||||||
}
|
|
||||||
|
|
||||||
essential_names curr_conc_names;
|
|
||||||
|
|
||||||
curr_conc_names[POET_SOL] = dump_solution_names(id);
|
|
||||||
init_values[id] = curr_conc_names;
|
|
||||||
|
|
||||||
// auto pair =
|
|
||||||
// std::make_pair(curr_conc_names,
|
|
||||||
// this->get_essential_values_init(id,
|
|
||||||
// curr_conc_names[0]));
|
|
||||||
// raw_initials[id] = pair;
|
|
||||||
}
|
|
||||||
|
|
||||||
std::vector<std::string> union_solution_names;
|
|
||||||
|
|
||||||
for (const auto &[id, val] : init_values) {
|
|
||||||
union_solution_names =
|
|
||||||
unionStringVectors(union_solution_names, val[POET_SOL]);
|
|
||||||
}
|
|
||||||
|
|
||||||
for (auto &[id, val] : init_values) {
|
|
||||||
dump_reactant_names(id, union_solution_names, val);
|
|
||||||
|
|
||||||
std::vector<std::string> solution_order(val[POET_SOL].begin() + 3,
|
|
||||||
val[POET_SOL].end());
|
|
||||||
|
|
||||||
auto pair = std::make_pair(
|
|
||||||
val, this->get_essential_values_init(id, solution_order));
|
|
||||||
raw_initials[id] = pair;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
std::string PhreeqcInit::subExchangeName(std::string name) {
|
|
||||||
for (const auto &species : this->PhreeqcPtr->Get_species_list()) {
|
|
||||||
const std::string &species_name = species.s->name;
|
|
||||||
// check if `name` is a prefix of `species_name`
|
|
||||||
if (species_name.compare(0, name.size(), name) == 0) {
|
|
||||||
return species_name;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
return name;
|
|
||||||
};
|
|
||||||
|
|
||||||
PhreeqcInit::essential_names PhreeqcInit::union_raws(const RawMap &raws) {
|
|
||||||
PhreeqcInit::essential_names names;
|
|
||||||
for (const auto &[sol_id, val] : raws) {
|
|
||||||
for (std::size_t i = 0; i < names.size(); i++) {
|
|
||||||
names[i] = unionStringVectors(names[i], val.first[i]);
|
|
||||||
}
|
|
||||||
|
|
||||||
for (auto &specie : names[POET_EXCH]) {
|
|
||||||
specie = this->subExchangeName(specie);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
return names;
|
|
||||||
}
|
|
||||||
|
|
||||||
std::vector<std::vector<double>>
|
|
||||||
PhreeqcInit::conc_from_essentials(const RawMap &raws,
|
|
||||||
const essential_names &names) {
|
|
||||||
std::vector<std::vector<double>> values;
|
|
||||||
|
|
||||||
for (const auto &[key, val] : raws) {
|
|
||||||
std::vector<double> combined_conc_vec;
|
|
||||||
|
|
||||||
for (std::size_t i = 0, offset = 0; i < names.size(); i++) {
|
|
||||||
std::vector<double> union_vec =
|
|
||||||
createConcVector(val.first[i], val.second.begin() + offset, names[i]);
|
|
||||||
|
|
||||||
combined_conc_vec.insert(combined_conc_vec.end(), union_vec.begin(),
|
|
||||||
union_vec.end());
|
|
||||||
|
|
||||||
offset += val.first[i].size();
|
|
||||||
}
|
|
||||||
|
|
||||||
values.push_back(combined_conc_vec);
|
|
||||||
}
|
|
||||||
|
|
||||||
return values;
|
|
||||||
}
|
|
||||||
|
|
||||||
PhreeqcInit::PhreeqcMat PhreeqcInit::buildPhreeqcMat() {
|
|
||||||
std::vector<std::vector<double>> values;
|
|
||||||
|
|
||||||
const PhreeqcInit::essential_names ess_names =
|
|
||||||
this->union_raws(this->raw_initials);
|
|
||||||
|
|
||||||
const std::vector<std::vector<double>> conc_values =
|
|
||||||
this->conc_from_essentials(this->raw_initials, ess_names);
|
|
||||||
|
|
||||||
std::vector<std::string> conc_names;
|
|
||||||
|
|
||||||
for (const auto &i : ess_names) {
|
|
||||||
conc_names.insert(conc_names.end(), i.begin(), i.end());
|
|
||||||
}
|
|
||||||
|
|
||||||
std::vector<int> ids(this->raw_initials.size());
|
|
||||||
|
|
||||||
std::transform(this->raw_initials.begin(), this->raw_initials.end(),
|
|
||||||
ids.begin(), [](const auto &pair) { return pair.first; });
|
|
||||||
|
|
||||||
return {.names = conc_names, .ids = ids, .values = conc_values};
|
|
||||||
}
|
|
||||||
|
|
||||||
PhreeqcInit::ModulesArray
|
|
||||||
PhreeqcInit::getModuleSizes(const std::vector<int> &cell_ids) {
|
|
||||||
ModulesArray module_sizes;
|
|
||||||
RawMap raws;
|
|
||||||
|
|
||||||
for (const auto &id : cell_ids) {
|
|
||||||
raws[id] = this->raw_initials[id];
|
|
||||||
}
|
|
||||||
|
|
||||||
const PhreeqcInit::essential_names ess_names = this->union_raws(raws);
|
|
||||||
|
|
||||||
std::transform(ess_names.begin(), ess_names.end(), module_sizes.begin(),
|
|
||||||
[](const auto &mod) { return mod.size(); });
|
|
||||||
|
|
||||||
return module_sizes;
|
|
||||||
}
|
|
||||||
|
|
||||||
// IPhreeqcPOET::PhreeqcMat
|
|
||||||
// IPhreeqcPOET::getPhreeqcMat(const std::vector<int> &use_ids) {
|
|
||||||
// std::vector<std::vector<double>> values;
|
|
||||||
// RawMap raws;
|
|
||||||
|
|
||||||
// for (const auto &id : use_ids) {
|
|
||||||
// raws[id] = this->raw_initials[id];
|
|
||||||
// }
|
|
||||||
|
|
||||||
// const IPhreeqcPOET::essential_names ess_names = this->union_raws(raws);
|
|
||||||
|
|
||||||
// const std::vector<std::vector<double>> conc_values =
|
|
||||||
// this->conc_from_essentials(raws, ess_names);
|
|
||||||
|
|
||||||
// std::vector<std::string> conc_names;
|
|
||||||
|
|
||||||
// for (const auto &i : ess_names) {
|
|
||||||
// conc_names.insert(conc_names.end(), i.begin(), i.end());
|
|
||||||
// }
|
|
||||||
|
|
||||||
// return {.names = conc_names, .ids = use_ids, .values = conc_values};
|
|
||||||
// }
|
|
||||||
@ -1,255 +0,0 @@
|
|||||||
#include "PhreeqcInit.hpp"
|
|
||||||
#include "Wrapper/EquilibriumWrapper.hpp"
|
|
||||||
#include "Wrapper/ExchangeWrapper.hpp"
|
|
||||||
#include "Wrapper/KineticWrapper.hpp"
|
|
||||||
#include "Wrapper/SolutionWrapper.hpp"
|
|
||||||
#include "Wrapper/SurfaceWrapper.hpp"
|
|
||||||
#include "global_structures.h"
|
|
||||||
|
|
||||||
#include <cassert>
|
|
||||||
#include <cstddef>
|
|
||||||
#include <set>
|
|
||||||
#include <span>
|
|
||||||
#include <string>
|
|
||||||
#include <vector>
|
|
||||||
|
|
||||||
std::vector<std::string> PhreeqcInit::dump_solution_names(int cell_number) {
|
|
||||||
constexpr std::size_t SKIP_H_O_CB = 3;
|
|
||||||
|
|
||||||
std::vector<std::string> placeholder;
|
|
||||||
|
|
||||||
return find_all_valence_states(
|
|
||||||
SolutionWrapper::names(this->Get_solution(cell_number), placeholder),
|
|
||||||
SKIP_H_O_CB);
|
|
||||||
}
|
|
||||||
|
|
||||||
void PhreeqcInit::dump_reactant_names(
|
|
||||||
int cell_number, const std::vector<std::string> union_sol_names,
|
|
||||||
essential_names &names) {
|
|
||||||
// Exchange
|
|
||||||
if (this->Get_exchange(cell_number) != NULL) {
|
|
||||||
auto *exchange = this->Get_exchange(cell_number);
|
|
||||||
|
|
||||||
std::vector<std::string> exch_formulas;
|
|
||||||
names[POET_EXCH] = ExchangeWrapper::names(exchange, exch_formulas);
|
|
||||||
|
|
||||||
this->exchanger[cell_number] = exch_formulas;
|
|
||||||
}
|
|
||||||
|
|
||||||
// Kinetics
|
|
||||||
if (this->Get_kinetic(cell_number) != NULL) {
|
|
||||||
auto *kinetic = this->Get_kinetic(cell_number);
|
|
||||||
|
|
||||||
std::vector<std::string> kin_comps;
|
|
||||||
|
|
||||||
names[POET_KIN] = KineticWrapper::names(kinetic, kin_comps);
|
|
||||||
|
|
||||||
this->kinetics[cell_number] = kin_comps;
|
|
||||||
}
|
|
||||||
|
|
||||||
// PPassemblage
|
|
||||||
if (this->Get_equilibrium(cell_number) != NULL) {
|
|
||||||
auto *equilibrium = this->Get_equilibrium(cell_number);
|
|
||||||
|
|
||||||
std::vector<std::string> equ_names;
|
|
||||||
|
|
||||||
names[POET_EQUIL] = EquilibriumWrapper::names(equilibrium, equ_names);
|
|
||||||
|
|
||||||
this->equilibrium[cell_number] = equ_names;
|
|
||||||
}
|
|
||||||
|
|
||||||
// Surface
|
|
||||||
if (this->Get_surface(cell_number) != NULL) {
|
|
||||||
auto *surface = this->Get_surface(cell_number);
|
|
||||||
|
|
||||||
if (this->surface_primaries.empty()) {
|
|
||||||
// this is fixed! Always add H and O
|
|
||||||
this->surface_primaries.insert("H");
|
|
||||||
this->surface_primaries.insert("O");
|
|
||||||
|
|
||||||
for (std::size_t i = 3; i < union_sol_names.size(); i++) {
|
|
||||||
const auto master_primary = this->PhreeqcPtr->master_bsearch_primary(
|
|
||||||
union_sol_names[i].c_str());
|
|
||||||
if (master_primary != NULL) {
|
|
||||||
this->surface_primaries.insert(master_primary->elt->name);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
std::vector<std::string> comp_formulas;
|
|
||||||
std::vector<std::string> charge_names;
|
|
||||||
|
|
||||||
names[POET_SURF] = SurfaceWrapper::names(surface, this->surface_primaries,
|
|
||||||
comp_formulas, charge_names);
|
|
||||||
|
|
||||||
this->surface_comps[cell_number] = comp_formulas;
|
|
||||||
this->surface_charge[cell_number] = charge_names;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
std::vector<std::string> PhreeqcInit::find_all_valence_states(
|
|
||||||
const std::vector<std::string> &solution_names, const std::size_t offset) {
|
|
||||||
std::vector<std::string> solution_with_valences(
|
|
||||||
solution_names.begin(), solution_names.begin() + offset);
|
|
||||||
|
|
||||||
// to avoid duplicates store already evaluated master species
|
|
||||||
std::set<std::string> master_species_found;
|
|
||||||
|
|
||||||
for (std::size_t i = offset; i < solution_names.size(); i++) {
|
|
||||||
const auto master_primary =
|
|
||||||
this->PhreeqcPtr->master_bsearch_primary(solution_names[i].c_str());
|
|
||||||
|
|
||||||
assert(master_primary != NULL);
|
|
||||||
|
|
||||||
const std::string master_species(master_primary->elt->name);
|
|
||||||
|
|
||||||
// in case we already found valences for this master species we skip it
|
|
||||||
if (master_species_found.contains(master_species)) {
|
|
||||||
continue;
|
|
||||||
}
|
|
||||||
|
|
||||||
master_species_found.insert(master_species);
|
|
||||||
|
|
||||||
bool has_valences = false;
|
|
||||||
std::size_t inner_loop_j = 0;
|
|
||||||
std::size_t last_valence = inner_loop_j;
|
|
||||||
|
|
||||||
// we HAVE to assume master species are already sorted!
|
|
||||||
// loop over all known master species
|
|
||||||
for (inner_loop_j = 0; inner_loop_j < this->PhreeqcPtr->master.size();
|
|
||||||
inner_loop_j++) {
|
|
||||||
std::string curr_master_species(
|
|
||||||
this->PhreeqcPtr->master[inner_loop_j]->elt->name);
|
|
||||||
|
|
||||||
if (curr_master_species == master_species) {
|
|
||||||
last_valence = inner_loop_j;
|
|
||||||
|
|
||||||
while (last_valence < this->PhreeqcPtr->master.size() - 1) {
|
|
||||||
std::string next_master_species(
|
|
||||||
this->PhreeqcPtr->master[last_valence + 1]->elt->name);
|
|
||||||
|
|
||||||
// check if the next name starts with the current master species
|
|
||||||
if (next_master_species.compare(0, master_species.size(),
|
|
||||||
master_species) != 0) {
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
|
|
||||||
// check if the next character is an opening parenthesis
|
|
||||||
if (next_master_species[master_species.size()] != '(') {
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
|
|
||||||
// if this is all true we have a valence state
|
|
||||||
has_valences = true;
|
|
||||||
last_valence++;
|
|
||||||
}
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// in case we found valences we add them to the solution
|
|
||||||
if (has_valences) {
|
|
||||||
// skip the master species and only add the valences
|
|
||||||
for (std::size_t k = inner_loop_j + 1; k <= last_valence; k++) {
|
|
||||||
solution_with_valences.push_back(
|
|
||||||
this->PhreeqcPtr->master[k]->elt->name);
|
|
||||||
}
|
|
||||||
continue;
|
|
||||||
}
|
|
||||||
|
|
||||||
// otherwise we just add the master species without any valences
|
|
||||||
solution_with_valences.push_back(master_species);
|
|
||||||
}
|
|
||||||
|
|
||||||
return solution_with_valences;
|
|
||||||
}
|
|
||||||
|
|
||||||
std::vector<double>
|
|
||||||
PhreeqcInit::get_essential_values_init(std::size_t cell_number,
|
|
||||||
const std::vector<std::string> &order) {
|
|
||||||
std::vector<double> essentials;
|
|
||||||
|
|
||||||
// Solutions
|
|
||||||
if (this->Get_solution(cell_number) != NULL) {
|
|
||||||
SolutionWrapper solution(this->Get_solution(cell_number), order);
|
|
||||||
std::vector<double> sol_values(solution.size());
|
|
||||||
std::span<double> sol_span{sol_values};
|
|
||||||
|
|
||||||
solution.get(sol_span);
|
|
||||||
|
|
||||||
essentials.insert(essentials.end(), sol_values.begin(), sol_values.end());
|
|
||||||
}
|
|
||||||
|
|
||||||
// Exchange
|
|
||||||
if (this->Get_exchange(cell_number) != NULL) {
|
|
||||||
ExchangeWrapper exchange(this->Get_exchange(cell_number),
|
|
||||||
this->exchanger[cell_number]);
|
|
||||||
|
|
||||||
std::vector<double> exch_values(exchange.size());
|
|
||||||
std::span<double> exch_span{exch_values};
|
|
||||||
|
|
||||||
exchange.get(exch_span);
|
|
||||||
|
|
||||||
essentials.insert(essentials.end(), exch_values.begin(), exch_values.end());
|
|
||||||
}
|
|
||||||
|
|
||||||
// Kinetics
|
|
||||||
if (this->Get_kinetic(cell_number) != NULL) {
|
|
||||||
KineticWrapper kinetic(this->Get_kinetic(cell_number),
|
|
||||||
this->kinetics[cell_number]);
|
|
||||||
|
|
||||||
std::vector<double> kin_values(kinetic.size());
|
|
||||||
std::span<double> kin_span{kin_values};
|
|
||||||
|
|
||||||
kinetic.get(kin_span);
|
|
||||||
|
|
||||||
essentials.insert(essentials.end(), kin_values.begin(), kin_values.end());
|
|
||||||
}
|
|
||||||
|
|
||||||
// PPassemblage
|
|
||||||
if (this->Get_equilibrium(cell_number) != NULL) {
|
|
||||||
EquilibriumWrapper equilibrium(this->Get_equilibrium(cell_number),
|
|
||||||
this->equilibrium[cell_number]);
|
|
||||||
|
|
||||||
std::vector<double> equ_values(equilibrium.size());
|
|
||||||
std::span<double> equ_span{equ_values};
|
|
||||||
|
|
||||||
equilibrium.get(equ_span);
|
|
||||||
|
|
||||||
essentials.insert(essentials.end(), equ_values.begin(), equ_values.end());
|
|
||||||
}
|
|
||||||
|
|
||||||
// Surface
|
|
||||||
if (this->Get_surface(cell_number) != NULL) {
|
|
||||||
|
|
||||||
SurfaceWrapper surface(
|
|
||||||
this->Get_surface(cell_number), this->surface_primaries,
|
|
||||||
this->surface_comps[cell_number], this->surface_charge[cell_number]);
|
|
||||||
|
|
||||||
std::vector<double> surf_values(surface.size());
|
|
||||||
std::span<double> surf_span{surf_values};
|
|
||||||
|
|
||||||
surface.get(surf_span);
|
|
||||||
|
|
||||||
essentials.insert(essentials.end(), surf_values.begin(), surf_values.end());
|
|
||||||
}
|
|
||||||
|
|
||||||
return essentials;
|
|
||||||
}
|
|
||||||
|
|
||||||
std::map<int, std::string> PhreeqcInit::raw_dumps() {
|
|
||||||
std::map<int, std::string> dumps;
|
|
||||||
|
|
||||||
this->SetDumpStringOn(true);
|
|
||||||
|
|
||||||
for (const auto &[sol_id, unused] : this->raw_initials) {
|
|
||||||
std::string call_string =
|
|
||||||
"DUMP\n -cells " + std::to_string(sol_id) + "\nEND\n";
|
|
||||||
this->RunString(call_string.c_str());
|
|
||||||
dumps[sol_id] = this->GetDumpString();
|
|
||||||
}
|
|
||||||
|
|
||||||
this->SetDumpStringOn(false);
|
|
||||||
|
|
||||||
return dumps;
|
|
||||||
}
|
|
||||||
247
poet/src/PhreeqcMatrix/Access.cpp
Normal file
247
poet/src/PhreeqcMatrix/Access.cpp
Normal file
@ -0,0 +1,247 @@
|
|||||||
|
#include <PhreeqcMatrix.hpp>
|
||||||
|
#include <algorithm>
|
||||||
|
#include <cstddef>
|
||||||
|
#include <limits>
|
||||||
|
#include <map>
|
||||||
|
#include <set>
|
||||||
|
#include <vector>
|
||||||
|
|
||||||
|
PhreeqcMatrix PhreeqcMatrix::subset(const std::vector<int> &indices) const {
|
||||||
|
PhreeqcMatrix result;
|
||||||
|
|
||||||
|
result._m_pqc = _m_pqc;
|
||||||
|
|
||||||
|
for (const auto &index : indices) {
|
||||||
|
result._m_map[index] = _m_map.at(index);
|
||||||
|
result._m_internal_names[index] = _m_internal_names.at(index);
|
||||||
|
}
|
||||||
|
|
||||||
|
result.remove_NaNs();
|
||||||
|
|
||||||
|
return result;
|
||||||
|
}
|
||||||
|
|
||||||
|
PhreeqcMatrix PhreeqcMatrix::erase(const std::vector<int> &indices) const {
|
||||||
|
PhreeqcMatrix result(*this);
|
||||||
|
|
||||||
|
for (const auto &index : indices) {
|
||||||
|
result._m_map.erase(index);
|
||||||
|
result._m_internal_names.erase(index);
|
||||||
|
}
|
||||||
|
|
||||||
|
result.remove_NaNs();
|
||||||
|
|
||||||
|
return result;
|
||||||
|
}
|
||||||
|
|
||||||
|
// PhreeqcMatrix &PhreeqcMatrix::subset(const std::vector<int> &indices) {
|
||||||
|
// std::map<int, std::vector<element>> new_map;
|
||||||
|
// std::map<int, std::vector<base_names>> new_internal_names;
|
||||||
|
|
||||||
|
// for (const auto &index : indices) {
|
||||||
|
// new_map[index] = _m_map.at(index);
|
||||||
|
// new_internal_names[index] = _m_internal_names.at(index);
|
||||||
|
// }
|
||||||
|
|
||||||
|
// this->_m_map = new_map;
|
||||||
|
// this->_m_internal_names = new_internal_names;
|
||||||
|
|
||||||
|
// this->remove_NaNs();
|
||||||
|
|
||||||
|
// return *this;
|
||||||
|
// }
|
||||||
|
|
||||||
|
// PhreeqcMatrix &PhreeqcMatrix::erase(const std::vector<int> &indices) {
|
||||||
|
// for (const auto &index : indices) {
|
||||||
|
// _m_map.erase(index);
|
||||||
|
// _m_internal_names.erase(index);
|
||||||
|
// }
|
||||||
|
|
||||||
|
// this->remove_NaNs();
|
||||||
|
|
||||||
|
// return *this;
|
||||||
|
// }
|
||||||
|
|
||||||
|
PhreeqcMatrix::STLExport PhreeqcMatrix::get(VectorExportType type,
|
||||||
|
bool include_id) const {
|
||||||
|
const auto &first_element = _m_map.begin()->second;
|
||||||
|
|
||||||
|
const std::size_t cols = first_element.size();
|
||||||
|
const std::size_t rows = _m_map.size();
|
||||||
|
const std::size_t total_elements = cols * rows;
|
||||||
|
|
||||||
|
STLExport result;
|
||||||
|
|
||||||
|
if (include_id) {
|
||||||
|
result.names.push_back("ID");
|
||||||
|
}
|
||||||
|
|
||||||
|
for (const auto &element : first_element) {
|
||||||
|
if (element.type != PhreeqcMatrix::PhreeqcComponent::SOLUTION) {
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
result.names.push_back(element.name);
|
||||||
|
}
|
||||||
|
|
||||||
|
for (std::size_t component = 1;
|
||||||
|
component <=
|
||||||
|
static_cast<std::size_t>(PhreeqcMatrix::PhreeqcComponent::SURFACE_COMPS);
|
||||||
|
component++) {
|
||||||
|
std::vector<std::string> values;
|
||||||
|
for (const auto &[id, elements] : _m_map) {
|
||||||
|
std::vector<std::string> names;
|
||||||
|
for (const auto &element : elements) {
|
||||||
|
if (static_cast<std::size_t>(element.type) == component) {
|
||||||
|
names.push_back(element.name);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
std::vector<std::string> union_names;
|
||||||
|
std::set_union(values.begin(), values.end(), names.begin(), names.end(),
|
||||||
|
std::back_inserter(union_names));
|
||||||
|
values = union_names;
|
||||||
|
}
|
||||||
|
result.names.insert(result.names.end(), values.begin(), values.end());
|
||||||
|
}
|
||||||
|
|
||||||
|
result.values.reserve(total_elements + (include_id ? rows : 0));
|
||||||
|
|
||||||
|
if (type == VectorExportType::COLUMN_MAJOR) {
|
||||||
|
std::size_t column_index = 0;
|
||||||
|
|
||||||
|
if (include_id) {
|
||||||
|
for (const auto &[id, _] : _m_map) {
|
||||||
|
result.values.push_back(id);
|
||||||
|
}
|
||||||
|
column_index++;
|
||||||
|
}
|
||||||
|
|
||||||
|
for (; column_index < result.names.size(); column_index++) {
|
||||||
|
for (const auto &[_, elements] : _m_map) {
|
||||||
|
double value_to_add = std::numeric_limits<double>::quiet_NaN();
|
||||||
|
for (const auto &curr_element : elements) {
|
||||||
|
const std::string &curr_element_name = curr_element.name;
|
||||||
|
|
||||||
|
if (curr_element_name == result.names[column_index]) {
|
||||||
|
value_to_add = curr_element.value;
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
result.values.push_back(value_to_add);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
return result;
|
||||||
|
}
|
||||||
|
|
||||||
|
for (const auto &[id, element_vec] : _m_map) {
|
||||||
|
std::size_t column_index = 0;
|
||||||
|
if (include_id) {
|
||||||
|
result.values.push_back(id);
|
||||||
|
column_index++;
|
||||||
|
}
|
||||||
|
|
||||||
|
for (; column_index < result.names.size(); column_index++) {
|
||||||
|
double value_to_add = std::numeric_limits<double>::quiet_NaN();
|
||||||
|
for (const auto &curr_element : element_vec) {
|
||||||
|
const std::string &curr_element_name = curr_element.name;
|
||||||
|
|
||||||
|
if (curr_element_name == result.names[column_index]) {
|
||||||
|
value_to_add = curr_element.value;
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
result.values.push_back(value_to_add);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
return result;
|
||||||
|
}
|
||||||
|
|
||||||
|
std::vector<std::string>
|
||||||
|
PhreeqcMatrix::getSolutionNames(bool include_h_o_charge) const {
|
||||||
|
std::vector<std::string> names;
|
||||||
|
|
||||||
|
if (include_h_o_charge) {
|
||||||
|
names.push_back("H");
|
||||||
|
names.push_back("O");
|
||||||
|
names.push_back("Charge");
|
||||||
|
}
|
||||||
|
|
||||||
|
const auto &first_element = _m_map.begin()->second;
|
||||||
|
|
||||||
|
for (std::size_t i = 3; i < _m_map.begin()->second.size(); i++) {
|
||||||
|
// assuming the element vector always starts with the solution components
|
||||||
|
if (first_element[i].type != PhreeqcComponent::SOLUTION) {
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
|
||||||
|
names.push_back(first_element[i].name);
|
||||||
|
}
|
||||||
|
|
||||||
|
return names;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <PhreeqcMatrix::base_names::Components comp>
|
||||||
|
static inline std::vector<std::string> get_names_from_internal_vector(
|
||||||
|
const std::map<int, std::vector<PhreeqcMatrix::base_names>> &map, int id) {
|
||||||
|
|
||||||
|
const auto &it = map.find(id);
|
||||||
|
|
||||||
|
if (it == map.end()) {
|
||||||
|
return {};
|
||||||
|
}
|
||||||
|
|
||||||
|
std::vector<std::string> names;
|
||||||
|
|
||||||
|
for (const auto &element : it->second) {
|
||||||
|
if (element.type == comp) {
|
||||||
|
names.push_back(element.name);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
return names;
|
||||||
|
}
|
||||||
|
|
||||||
|
std::vector<std::string> PhreeqcMatrix::getExchanger(int id) const {
|
||||||
|
return get_names_from_internal_vector<base_names::Components::EXCHANGER>(
|
||||||
|
this->_m_internal_names, id);
|
||||||
|
}
|
||||||
|
|
||||||
|
std::vector<std::string> PhreeqcMatrix::getKineticsNames(int id) const {
|
||||||
|
return get_names_from_internal_vector<base_names::Components::KINETICS>(
|
||||||
|
this->_m_internal_names, id);
|
||||||
|
}
|
||||||
|
|
||||||
|
std::vector<std::string> PhreeqcMatrix::getEquilibriumNames(int id) const {
|
||||||
|
return get_names_from_internal_vector<base_names::Components::EQUILIBRIUM>(
|
||||||
|
this->_m_internal_names, id);
|
||||||
|
}
|
||||||
|
|
||||||
|
std::vector<std::string> PhreeqcMatrix::getSurfaceCompNames(int id) const {
|
||||||
|
return get_names_from_internal_vector<base_names::Components::SURACE_COMP>(
|
||||||
|
this->_m_internal_names, id);
|
||||||
|
}
|
||||||
|
|
||||||
|
std::vector<std::string> PhreeqcMatrix::getSurfaceChargeNames(int id) const {
|
||||||
|
return get_names_from_internal_vector<base_names::Components::SURFACE_CHARGE>(
|
||||||
|
this->_m_internal_names, id);
|
||||||
|
}
|
||||||
|
|
||||||
|
std::vector<std::string> PhreeqcMatrix::getSolutionPrimaries() const {
|
||||||
|
return std::vector<std::string>(_m_surface_primaries.begin(),
|
||||||
|
_m_surface_primaries.end());
|
||||||
|
}
|
||||||
|
|
||||||
|
double PhreeqcMatrix::operator()(int cell_id, const std::string &name) const {
|
||||||
|
const auto &elements = _m_map.at(cell_id);
|
||||||
|
|
||||||
|
const auto it = std::find_if(
|
||||||
|
elements.begin(), elements.end(),
|
||||||
|
[&name](const element &element) { return element.name == name; });
|
||||||
|
|
||||||
|
if (it == elements.end()) {
|
||||||
|
throw std::runtime_error("Element not found");
|
||||||
|
}
|
||||||
|
|
||||||
|
return it->value;
|
||||||
|
}
|
||||||
48
poet/src/PhreeqcMatrix/Ctor.cpp
Normal file
48
poet/src/PhreeqcMatrix/Ctor.cpp
Normal file
@ -0,0 +1,48 @@
|
|||||||
|
#include "IPhreeqc.hpp"
|
||||||
|
#include "PhreeqcMatrix.hpp"
|
||||||
|
|
||||||
|
#include <Phreeqc.h>
|
||||||
|
#include <Solution.h>
|
||||||
|
#include <cmath>
|
||||||
|
#include <string>
|
||||||
|
|
||||||
|
PhreeqcMatrix::PhreeqcMatrix(const std::string &database,
|
||||||
|
const std::string &input_script)
|
||||||
|
: _m_database(database) {
|
||||||
|
this->_m_pqc = std::make_shared<IPhreeqc>();
|
||||||
|
|
||||||
|
this->_m_pqc->LoadDatabaseString(database.c_str());
|
||||||
|
this->_m_pqc->RunString(input_script.c_str());
|
||||||
|
|
||||||
|
this->initialize();
|
||||||
|
}
|
||||||
|
|
||||||
|
PhreeqcMatrix::PhreeqcMatrix(const PhreeqcMatrix &other)
|
||||||
|
: _m_map(other._m_map), _m_internal_names(other._m_internal_names),
|
||||||
|
_m_surface_primaries(other._m_surface_primaries), _m_pqc(other._m_pqc),
|
||||||
|
_m_database(other._m_database) {}
|
||||||
|
|
||||||
|
PhreeqcMatrix::PhreeqcMatrix(PhreeqcMatrix &&other)
|
||||||
|
: _m_map(other._m_map), _m_internal_names(other._m_internal_names),
|
||||||
|
_m_surface_primaries(other._m_surface_primaries), _m_pqc(other._m_pqc),
|
||||||
|
_m_database(other._m_database) {}
|
||||||
|
|
||||||
|
PhreeqcMatrix &PhreeqcMatrix::operator=(const PhreeqcMatrix &other) {
|
||||||
|
_m_map = other._m_map;
|
||||||
|
_m_internal_names = other._m_internal_names;
|
||||||
|
_m_surface_primaries = other._m_surface_primaries;
|
||||||
|
_m_pqc = other._m_pqc;
|
||||||
|
_m_database = other._m_database;
|
||||||
|
|
||||||
|
return *this;
|
||||||
|
}
|
||||||
|
|
||||||
|
PhreeqcMatrix &PhreeqcMatrix::operator=(PhreeqcMatrix &&other) {
|
||||||
|
_m_map = other._m_map;
|
||||||
|
_m_internal_names = other._m_internal_names;
|
||||||
|
_m_surface_primaries = other._m_surface_primaries;
|
||||||
|
_m_pqc = other._m_pqc;
|
||||||
|
_m_database = other._m_database;
|
||||||
|
|
||||||
|
return *this;
|
||||||
|
}
|
||||||
254
poet/src/PhreeqcMatrix/Init.cpp
Normal file
254
poet/src/PhreeqcMatrix/Init.cpp
Normal file
@ -0,0 +1,254 @@
|
|||||||
|
#include "PhreeqcMatrix.hpp"
|
||||||
|
|
||||||
|
#include "../Wrapper/EquilibriumWrapper.hpp"
|
||||||
|
#include "../Wrapper/ExchangeWrapper.hpp"
|
||||||
|
#include "../Wrapper/KineticWrapper.hpp"
|
||||||
|
#include "../Wrapper/SolutionWrapper.hpp"
|
||||||
|
#include "../Wrapper/SurfaceWrapper.hpp"
|
||||||
|
|
||||||
|
#include <Phreeqc.h>
|
||||||
|
#include <Solution.h>
|
||||||
|
#include <cmath>
|
||||||
|
#include <cstddef>
|
||||||
|
#include <iterator>
|
||||||
|
#include <map>
|
||||||
|
#include <set>
|
||||||
|
#include <string>
|
||||||
|
#include <utility>
|
||||||
|
#include <vector>
|
||||||
|
|
||||||
|
constexpr std::size_t SKIP_H_O_CB = 3;
|
||||||
|
|
||||||
|
static std::vector<std::string> dump_solution_names(cxxSolution *solution,
|
||||||
|
Phreeqc *phreeqc) {
|
||||||
|
std::vector<std::string> placeholder;
|
||||||
|
|
||||||
|
return phreeqc->find_all_valence_states(
|
||||||
|
SolutionWrapper::names(solution, placeholder), SKIP_H_O_CB);
|
||||||
|
}
|
||||||
|
|
||||||
|
template <enum PhreeqcMatrix::PhreeqcComponent comp, class T>
|
||||||
|
static void
|
||||||
|
base_add_to_element_vector(const T &wrapper,
|
||||||
|
const std::vector<std::string> &names,
|
||||||
|
std::vector<PhreeqcMatrix::element> &elements) {
|
||||||
|
std::vector<double> values(wrapper.size());
|
||||||
|
std::span<double> values_span(values);
|
||||||
|
|
||||||
|
wrapper.get(values_span);
|
||||||
|
|
||||||
|
for (std::size_t i = 0; i < names.size(); ++i) {
|
||||||
|
elements.push_back({names[i], comp, values[i]});
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
template <enum PhreeqcMatrix::PhreeqcComponent comp, class Wrapper_t, class T>
|
||||||
|
static void ex_ki_eq_add_to_element_vector(
|
||||||
|
T *phreeqc_component, std::vector<PhreeqcMatrix::element> &elements,
|
||||||
|
std::vector<PhreeqcMatrix::base_names> &b_names) {
|
||||||
|
if (phreeqc_component == NULL) {
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
std::vector<std::string> ctor_names;
|
||||||
|
std::vector<std::string> names =
|
||||||
|
Wrapper_t::names(phreeqc_component, ctor_names);
|
||||||
|
|
||||||
|
for (const auto &name : ctor_names) {
|
||||||
|
b_names.push_back(
|
||||||
|
{static_cast<PhreeqcMatrix::base_names::Components>(comp), name});
|
||||||
|
}
|
||||||
|
|
||||||
|
Wrapper_t wrapper(phreeqc_component, ctor_names);
|
||||||
|
|
||||||
|
base_add_to_element_vector<comp>(wrapper, names, elements);
|
||||||
|
}
|
||||||
|
|
||||||
|
static void
|
||||||
|
surface_add_to_element_vector(Phreeqc *phreeqc, cxxSurface *surface,
|
||||||
|
const std::vector<std::string> &solution_names,
|
||||||
|
std::vector<PhreeqcMatrix::element> &elements,
|
||||||
|
std::vector<PhreeqcMatrix::base_names> &b_names,
|
||||||
|
std::set<std::string> &surface_primaries) {
|
||||||
|
|
||||||
|
if (surface == NULL) {
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
// H and O are fixed surface primaries
|
||||||
|
surface_primaries.insert("H");
|
||||||
|
surface_primaries.insert("O");
|
||||||
|
|
||||||
|
for (std::size_t i = 3; i < solution_names.size(); i++) {
|
||||||
|
const auto master_primary =
|
||||||
|
phreeqc->master_bsearch_primary(solution_names[i].c_str());
|
||||||
|
if (master_primary != NULL) {
|
||||||
|
surface_primaries.insert(master_primary->elt->name);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
std::vector<std::string> comp_formulas;
|
||||||
|
std::vector<std::string> charge_names;
|
||||||
|
|
||||||
|
std::vector<std::string> surf_names = SurfaceWrapper::names(
|
||||||
|
surface, surface_primaries, comp_formulas, charge_names);
|
||||||
|
|
||||||
|
for (const auto &comp : comp_formulas) {
|
||||||
|
b_names.push_back({static_cast<PhreeqcMatrix::base_names::Components>(
|
||||||
|
PhreeqcMatrix::base_names::Components::SURACE_COMP),
|
||||||
|
comp});
|
||||||
|
}
|
||||||
|
|
||||||
|
for (const auto &charge : charge_names) {
|
||||||
|
b_names.push_back(
|
||||||
|
{static_cast<PhreeqcMatrix::base_names::Components>(
|
||||||
|
PhreeqcMatrix::base_names::Components::SURFACE_CHARGE),
|
||||||
|
charge});
|
||||||
|
}
|
||||||
|
|
||||||
|
SurfaceWrapper wrapper(surface, surface_primaries, comp_formulas,
|
||||||
|
charge_names);
|
||||||
|
|
||||||
|
base_add_to_element_vector<PhreeqcMatrix::PhreeqcComponent::SURFACE_COMPS>(
|
||||||
|
wrapper, surf_names, elements);
|
||||||
|
}
|
||||||
|
|
||||||
|
static std::pair<std::vector<PhreeqcMatrix::element>,
|
||||||
|
std::vector<PhreeqcMatrix::base_names>>
|
||||||
|
create_vector_from_phreeqc(Phreeqc *phreeqc, int id,
|
||||||
|
const std::vector<std::string> &solution_names,
|
||||||
|
std::set<std::string> &surface_primaries) {
|
||||||
|
std::vector<PhreeqcMatrix::element> elements;
|
||||||
|
std::vector<PhreeqcMatrix::base_names> b_names;
|
||||||
|
|
||||||
|
// Solution
|
||||||
|
SolutionWrapper sol_wrapper(
|
||||||
|
Utilities::Rxn_find(phreeqc->Get_Rxn_solution_map(), id),
|
||||||
|
{solution_names.begin() + SKIP_H_O_CB, solution_names.end()});
|
||||||
|
|
||||||
|
base_add_to_element_vector<PhreeqcMatrix::PhreeqcComponent::SOLUTION>(
|
||||||
|
sol_wrapper, solution_names, elements);
|
||||||
|
|
||||||
|
// keep track if exchange is present
|
||||||
|
const size_t vec_size_before_exchange = elements.size();
|
||||||
|
|
||||||
|
// Exchange
|
||||||
|
ex_ki_eq_add_to_element_vector<PhreeqcMatrix::PhreeqcComponent::EXCHANGE,
|
||||||
|
ExchangeWrapper>(
|
||||||
|
Utilities::Rxn_find(phreeqc->Get_Rxn_exchange_map(), id), elements,
|
||||||
|
b_names);
|
||||||
|
|
||||||
|
const bool has_exchange = elements.size() > vec_size_before_exchange;
|
||||||
|
|
||||||
|
// Kinetic
|
||||||
|
ex_ki_eq_add_to_element_vector<PhreeqcMatrix::PhreeqcComponent::KINETIC,
|
||||||
|
KineticWrapper>(
|
||||||
|
Utilities::Rxn_find(phreeqc->Get_Rxn_kinetics_map(), id), elements,
|
||||||
|
b_names);
|
||||||
|
// Equilibrium
|
||||||
|
ex_ki_eq_add_to_element_vector<PhreeqcMatrix::PhreeqcComponent::EQUILIBRIUM,
|
||||||
|
EquilibriumWrapper>(
|
||||||
|
Utilities::Rxn_find(phreeqc->Get_Rxn_pp_assemblage_map(), id), elements,
|
||||||
|
b_names);
|
||||||
|
|
||||||
|
// Surface
|
||||||
|
surface_add_to_element_vector(
|
||||||
|
phreeqc, Utilities::Rxn_find(phreeqc->Get_Rxn_surface_map(), id),
|
||||||
|
solution_names, elements, b_names, surface_primaries);
|
||||||
|
|
||||||
|
// substitute exchange names
|
||||||
|
if (has_exchange) {
|
||||||
|
for (auto ¤tElement : elements) {
|
||||||
|
if (currentElement.type == PhreeqcMatrix::PhreeqcComponent::EXCHANGE) {
|
||||||
|
for (const auto &species : phreeqc->Get_species_list()) {
|
||||||
|
const std::string &species_name = species.s->name;
|
||||||
|
// check if `name` is a prefix of `species_name`
|
||||||
|
if (species_name.compare(0, currentElement.name.size(),
|
||||||
|
currentElement.name) == 0) {
|
||||||
|
currentElement.name = species_name;
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
return {elements, b_names};
|
||||||
|
}
|
||||||
|
|
||||||
|
static std::vector<std::string> find_all_solutions(Phreeqc *phreeqc) {
|
||||||
|
std::vector<std::vector<std::string>> all_names;
|
||||||
|
|
||||||
|
for (auto &[id, solution] : phreeqc->Get_Rxn_solution_map()) {
|
||||||
|
if (id < 0) {
|
||||||
|
continue;
|
||||||
|
}
|
||||||
|
std::vector<std::string> names = dump_solution_names(&solution, phreeqc);
|
||||||
|
all_names.push_back(names);
|
||||||
|
}
|
||||||
|
|
||||||
|
std::vector<std::string> union_names;
|
||||||
|
|
||||||
|
for (const auto &vec : all_names) {
|
||||||
|
std::vector<std::string> result;
|
||||||
|
std::set_union(union_names.begin(), union_names.end(), vec.begin(),
|
||||||
|
vec.end(), std::back_inserter(result));
|
||||||
|
|
||||||
|
union_names = result;
|
||||||
|
}
|
||||||
|
|
||||||
|
return union_names;
|
||||||
|
}
|
||||||
|
|
||||||
|
void PhreeqcMatrix::initialize() {
|
||||||
|
|
||||||
|
Phreeqc *phreeqc = this->_m_pqc->GetPhreeqcPtr();
|
||||||
|
|
||||||
|
if (!phreeqc->Get_Rxn_surface_map().empty()) {
|
||||||
|
this->_m_pqc->RunString("RUN_CELLS\n-cells 1\nEND");
|
||||||
|
}
|
||||||
|
|
||||||
|
std::vector<std::string> solutions = find_all_solutions(phreeqc);
|
||||||
|
|
||||||
|
for (auto &[id, solution] : phreeqc->Get_Rxn_solution_map()) {
|
||||||
|
if (id < 0) {
|
||||||
|
continue;
|
||||||
|
}
|
||||||
|
const auto &[elements, base_names] = create_vector_from_phreeqc(
|
||||||
|
phreeqc, id, solutions, this->_m_surface_primaries);
|
||||||
|
|
||||||
|
_m_map[id] = elements;
|
||||||
|
_m_internal_names[id] = base_names;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
void PhreeqcMatrix::remove_NaNs() {
|
||||||
|
std::vector<std::vector<std::string>> elements_to_remove;
|
||||||
|
|
||||||
|
for (const auto &[id, elements] : _m_map) {
|
||||||
|
elements_to_remove.push_back({});
|
||||||
|
std::vector<std::string> &curr_to_remove = elements_to_remove.back();
|
||||||
|
for (const auto &element : elements) {
|
||||||
|
if (std::isnan(element.value)) {
|
||||||
|
curr_to_remove.push_back(element.name);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// find intersection of all elements to remove
|
||||||
|
std::vector<std::string> intersection = elements_to_remove[0];
|
||||||
|
|
||||||
|
for (const auto &vec : elements_to_remove) {
|
||||||
|
std::set_intersection(intersection.begin(), intersection.end(), vec.begin(),
|
||||||
|
vec.end(), std::back_inserter(intersection));
|
||||||
|
}
|
||||||
|
|
||||||
|
for (auto &[id, elements] : _m_map) {
|
||||||
|
for (const auto &to_remove_element : intersection) {
|
||||||
|
elements.erase(std::remove_if(elements.begin(), elements.end(),
|
||||||
|
[to_remove_element](const element &el) {
|
||||||
|
return el.name == to_remove_element;
|
||||||
|
}),
|
||||||
|
elements.end());
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
53
poet/src/PhreeqcMatrix/Misc.cpp
Normal file
53
poet/src/PhreeqcMatrix/Misc.cpp
Normal file
@ -0,0 +1,53 @@
|
|||||||
|
#include "PhreeqcMatrix.hpp"
|
||||||
|
#include <cstddef>
|
||||||
|
|
||||||
|
std::vector<int> PhreeqcMatrix::getIds() const {
|
||||||
|
std::vector<int> ids;
|
||||||
|
for (const auto &[id, _] : _m_map) {
|
||||||
|
ids.push_back(id);
|
||||||
|
}
|
||||||
|
|
||||||
|
return ids;
|
||||||
|
}
|
||||||
|
|
||||||
|
// std::array<std::size_t, 5> PhreeqcMatrix::getComponentCount(int cell_id)
|
||||||
|
// const {
|
||||||
|
// std::array<std::size_t, 5> counts = {0, 0, 0, 0, 0};
|
||||||
|
|
||||||
|
// if (!this->_m_internal_names.contains(cell_id)) {
|
||||||
|
// return counts;
|
||||||
|
// }
|
||||||
|
|
||||||
|
// for (const auto &element : this->_m_internal_names.at(cell_id)) {
|
||||||
|
// counts[static_cast<std::size_t>(element.type)]++;
|
||||||
|
// }
|
||||||
|
|
||||||
|
// return counts;
|
||||||
|
// }
|
||||||
|
|
||||||
|
std::map<int, std::string> PhreeqcMatrix::getDumpStringsPQI() const {
|
||||||
|
std::map<int, std::string> dumps;
|
||||||
|
|
||||||
|
for (const auto &cell_id : this->getIds()) {
|
||||||
|
dumps[cell_id] = this->getDumpStringsPQI(cell_id);
|
||||||
|
}
|
||||||
|
|
||||||
|
return dumps;
|
||||||
|
}
|
||||||
|
|
||||||
|
std::string PhreeqcMatrix::getDumpStringsPQI(int cell_id) const {
|
||||||
|
this->_m_pqc->SetDumpStringOn(true);
|
||||||
|
|
||||||
|
const std::string call_string =
|
||||||
|
"DUMP\n -cells " + std::to_string(cell_id) + "\nEND\n";
|
||||||
|
|
||||||
|
this->_m_pqc->RunString(call_string.c_str());
|
||||||
|
|
||||||
|
const std::string dump_string = this->_m_pqc->GetDumpString();
|
||||||
|
|
||||||
|
this->_m_pqc->SetDumpStringOn(false);
|
||||||
|
|
||||||
|
return dump_string;
|
||||||
|
}
|
||||||
|
|
||||||
|
std::string PhreeqcMatrix::getDatabase() const { return _m_database; }
|
||||||
@ -23,7 +23,7 @@ std::vector<std::string> EquilibriumWrapper::EquilibriumCompWrapper::names(
|
|||||||
std::vector<std::string> names;
|
std::vector<std::string> names;
|
||||||
|
|
||||||
const std::string &comp_name = comp.Get_name();
|
const std::string &comp_name = comp.Get_name();
|
||||||
names.push_back(comp_name);
|
names.push_back(comp_name + "_eq");
|
||||||
names.push_back(comp_name + "_si");
|
names.push_back(comp_name + "_si");
|
||||||
|
|
||||||
return names;
|
return names;
|
||||||
|
|||||||
24
poet/test/CMakeLists.txt
Normal file
24
poet/test/CMakeLists.txt
Normal file
@ -0,0 +1,24 @@
|
|||||||
|
enable_testing()
|
||||||
|
|
||||||
|
add_executable(
|
||||||
|
poet_test
|
||||||
|
testPhreeqcMatrix.cpp
|
||||||
|
)
|
||||||
|
|
||||||
|
target_link_libraries(
|
||||||
|
poet_test
|
||||||
|
IPhreeqcPOET
|
||||||
|
GTest::gtest_main
|
||||||
|
)
|
||||||
|
|
||||||
|
target_include_directories(poet_test PRIVATE ${CMAKE_CURRENT_BINARY_DIR})
|
||||||
|
|
||||||
|
# read file and store in variable
|
||||||
|
file(READ "${PROJECT_SOURCE_DIR}/database/phreeqc.dat" POET_PHREEQCDAT_DB)
|
||||||
|
file(READ "${CMAKE_CURRENT_SOURCE_DIR}/barite_db.dat" POET_BARITE_DB)
|
||||||
|
file(READ "${CMAKE_CURRENT_SOURCE_DIR}/barite_het.pqi" POET_BARITE_PQI)
|
||||||
|
|
||||||
|
configure_file("${CMAKE_CURRENT_SOURCE_DIR}/testPhreeqcMatrix.hpp.in" "${CMAKE_CURRENT_BINARY_DIR}/testPhreeqcMatrix.hpp")
|
||||||
|
|
||||||
|
include(GoogleTest)
|
||||||
|
gtest_discover_tests(poet_test)
|
||||||
195
poet/test/barite_db.dat
Normal file
195
poet/test/barite_db.dat
Normal file
@ -0,0 +1,195 @@
|
|||||||
|
DATABASE
|
||||||
|
###########################
|
||||||
|
SOLUTION_MASTER_SPECIES
|
||||||
|
H H+ -1 H 1.008 # phreeqc/
|
||||||
|
H(0) H2 0 H # phreeqc/
|
||||||
|
H(1) H+ -1 0.0 # phreeqc/
|
||||||
|
E e- 0 0.0 0.0 # phreeqc/
|
||||||
|
O H2O 0 O 16.0 # phreeqc/
|
||||||
|
O(0) O2 0 O # phreeqc/
|
||||||
|
O(-2) H2O 0 0.0 # phreeqc/
|
||||||
|
Na Na+ 0 Na 22.9898 # phreeqc/
|
||||||
|
Ba Ba+2 0 Ba 137.34 # phreeqc/
|
||||||
|
Sr Sr+2 0 Sr 87.62 # phreeqc/
|
||||||
|
Cl Cl- 0 Cl 35.453 # phreeqc/
|
||||||
|
S SO4-2 0 SO4 32.064 # phreeqc/
|
||||||
|
S(6) SO4-2 0 SO4 # phreeqc/
|
||||||
|
S(-2) HS- 1 S # phreeqc/
|
||||||
|
SOLUTION_SPECIES
|
||||||
|
H+ = H+
|
||||||
|
-gamma 9 0
|
||||||
|
-dw 9.31e-09
|
||||||
|
# source: phreeqc
|
||||||
|
e- = e-
|
||||||
|
# source: phreeqc
|
||||||
|
H2O = H2O
|
||||||
|
# source: phreeqc
|
||||||
|
Na+ = Na+
|
||||||
|
-gamma 4.08 0.082
|
||||||
|
-dw 1.33e-09
|
||||||
|
-Vm 2.28 -4.38 -4.1 -0.586 0.09 4 0.3 52 -0.00333 0.566
|
||||||
|
# source: phreeqc
|
||||||
|
Ba+2 = Ba+2
|
||||||
|
-gamma 4 0.153
|
||||||
|
-dw 8.48e-10
|
||||||
|
-Vm 2.063 -10.06 1.9534 -2.36 0.4218 5 1.58 -12.03 -0.00835 1
|
||||||
|
# source: phreeqc
|
||||||
|
Sr+2 = Sr+2
|
||||||
|
-gamma 5.26 0.121
|
||||||
|
-dw 7.94e-10
|
||||||
|
-Vm -0.0157 -10.15 10.18 -2.36 0.86 5.26 0.859 -27 -0.0041 1.97
|
||||||
|
# source: phreeqc
|
||||||
|
Cl- = Cl-
|
||||||
|
-gamma 3.63 0.017
|
||||||
|
-dw 2.03e-09
|
||||||
|
-Vm 4.465 4.801 4.325 -2.847 1.748 0 -0.331 20.16 0 1
|
||||||
|
# source: phreeqc
|
||||||
|
SO4-2 = SO4-2
|
||||||
|
-gamma 5 -0.04
|
||||||
|
-dw 1.07e-09
|
||||||
|
-Vm 8 2.3 -46.04 6.245 3.82 0 0 0 0 1
|
||||||
|
# source: phreeqc
|
||||||
|
H2O = OH- + H+
|
||||||
|
-analytical_expression 293.29227 0.1360833 -10576.913 -123.73158 0 -6.996455e-05
|
||||||
|
-gamma 3.5 0
|
||||||
|
-dw 5.27e-09
|
||||||
|
-Vm -9.66 28.5 80 -22.9 1.89 0 1.09 0 0 1
|
||||||
|
# source: phreeqc
|
||||||
|
2 H2O = O2 + 4 H+ + 4 e-
|
||||||
|
-log_k -86.08
|
||||||
|
-delta_h 134.79 kcal
|
||||||
|
-dw 2.35e-09
|
||||||
|
-Vm 5.7889 6.3536 3.2528 -3.0417 -0.3943
|
||||||
|
# source: phreeqc
|
||||||
|
2 H+ + 2 e- = H2
|
||||||
|
-log_k -3.15
|
||||||
|
-delta_h -1.759 kcal
|
||||||
|
-dw 5.13e-09
|
||||||
|
-Vm 6.52 0.78 0.12
|
||||||
|
# source: phreeqc
|
||||||
|
SO4-2 + H+ = HSO4-
|
||||||
|
-log_k 1.988
|
||||||
|
-delta_h 3.85 kcal
|
||||||
|
-analytical_expression -56.889 0.006473 2307.9 19.8858
|
||||||
|
-dw 1.33e-09
|
||||||
|
-Vm 8.2 9.259 2.1108 -3.1618 1.1748 0 -0.3 15 0 1
|
||||||
|
# source: phreeqc
|
||||||
|
HS- = S-2 + H+
|
||||||
|
-log_k -12.918
|
||||||
|
-delta_h 12.1 kcal
|
||||||
|
-gamma 5 0
|
||||||
|
-dw 7.31e-10
|
||||||
|
# source: phreeqc
|
||||||
|
SO4-2 + 9 H+ + 8 e- = HS- + 4 H2O
|
||||||
|
-log_k 33.65
|
||||||
|
-delta_h -60.140 kcal
|
||||||
|
-gamma 3.5 0
|
||||||
|
-dw 1.73e-09
|
||||||
|
-Vm 5.0119 4.9799 3.4765 -2.9849 1.441
|
||||||
|
# source: phreeqc
|
||||||
|
HS- + H+ = H2S
|
||||||
|
-log_k 6.994
|
||||||
|
-delta_h -5.30 kcal
|
||||||
|
-analytical_expression -11.17 0.02386 3279
|
||||||
|
-dw 2.1e-09
|
||||||
|
-Vm 7.81 2.96 -0.46
|
||||||
|
# source: phreeqc
|
||||||
|
Na+ + OH- = NaOH
|
||||||
|
-log_k -10
|
||||||
|
# source: phreeqc
|
||||||
|
Na+ + SO4-2 = NaSO4-
|
||||||
|
-log_k 0.7
|
||||||
|
-delta_h 1.120 kcal
|
||||||
|
-gamma 5.4 0
|
||||||
|
-dw 1.33e-09
|
||||||
|
-Vm 1e-05 16.4 -0.0678 -1.05 4.14 0 6.86 0 0.0242 0.53
|
||||||
|
# source: phreeqc
|
||||||
|
Ba+2 + H2O = BaOH+ + H+
|
||||||
|
-log_k -13.47
|
||||||
|
-gamma 5 0
|
||||||
|
# source: phreeqc
|
||||||
|
Ba+2 + SO4-2 = BaSO4
|
||||||
|
-log_k 2.7
|
||||||
|
# source: phreeqc
|
||||||
|
Sr+2 + H2O = SrOH+ + H+
|
||||||
|
-log_k -13.29
|
||||||
|
-gamma 5 0
|
||||||
|
# source: phreeqc
|
||||||
|
Sr+2 + SO4-2 = SrSO4
|
||||||
|
-log_k 2.29
|
||||||
|
-delta_h 2.08 kcal
|
||||||
|
-Vm 6.791 -0.9666 6.13 -2.739 -0.001
|
||||||
|
# source: phreeqc
|
||||||
|
PHASES
|
||||||
|
Barite
|
||||||
|
BaSO4 = Ba+2 + SO4-2
|
||||||
|
-log_k -9.97
|
||||||
|
-delta_h 6.35 kcal
|
||||||
|
-analytical_expression -282.43 -0.08972 5822 113.08
|
||||||
|
-Vm 52.9
|
||||||
|
# source: phreeqc
|
||||||
|
# comment:
|
||||||
|
Celestite
|
||||||
|
SrSO4 = Sr+2 + SO4-2
|
||||||
|
-log_k -6.63
|
||||||
|
-delta_h -4.037 kcal
|
||||||
|
-analytical_expression -7.14 0.00611 75 0 0 -1.79e-05
|
||||||
|
-Vm 46.4
|
||||||
|
# source: phreeqc
|
||||||
|
# comment:
|
||||||
|
RATES
|
||||||
|
Celestite # Palandri & Kharaka 2004<--------------------------------change me
|
||||||
|
# PARM(1): reactive surface area
|
||||||
|
# am: acid mechanism, nm: neutral mechanism, bm: base mechanism
|
||||||
|
-start
|
||||||
|
10 sr_i = SR("Celestite") # saturation ratio, (-)<----------change me
|
||||||
|
20 moles = 0 # init target variable, (mol)
|
||||||
|
30 IF ((M <= 0) AND (sr_i < 1)) OR (sr_i = 1.0) THEN GOTO 310
|
||||||
|
40 sa = PARM(1) # reactive surface area, (m2)
|
||||||
|
|
||||||
|
100 r = 8.314462 # gas constant, (J K-1 mol-1)
|
||||||
|
110 dTi = (1 / TK) - (1 / 298.15) # (K-1)
|
||||||
|
120 ea_am = 23800 # activation energy am, (J mol-1)<-----------change me
|
||||||
|
130 ea_nm = 0 # activation energy nm, (J mol-1)<-----------change me
|
||||||
|
140 ea_bm = 0 # activation energy bm, (J mol-1)<-----------change me
|
||||||
|
150 log_k_am = -5.66 # reaction constant am<-------------------change me
|
||||||
|
rem log_k_nm = -99 # reaction constant nm<-------------------change me
|
||||||
|
rem log_k_bm = -99 # reaction constant bm<-------------------change me
|
||||||
|
180 n_am = 0.109 # H+ reaction order am<-----------------------change me
|
||||||
|
rem n_bm = 0 # H+ reaction order bm<-----------------------change me
|
||||||
|
200 am = (10 ^ log_k_am) * EXP(-ea_am * dTi / r) * ACT("H+") ^ n_am
|
||||||
|
rem nm = (10 ^ log_k_nm) * EXP(-ea_nm * dTi / r)
|
||||||
|
rem bm = (10 ^ log_k_bm) * EXP(-ea_bm * dTi / r) * ACT("H+") ^ n_bm
|
||||||
|
|
||||||
|
300 moles = sa * (am) * (1 - sr_i)
|
||||||
|
310 save moles * time
|
||||||
|
-end
|
||||||
|
|
||||||
|
Barite # Palandri & Kharaka 2004<-----------------------------------change me
|
||||||
|
# PARM(1): reactive surface area
|
||||||
|
# am: acid mechanism, nm: neutral mechanism, bm: base mechanism
|
||||||
|
-start
|
||||||
|
10 sr_i = SR("Barite") # saturation ratio, (-)<----------change me
|
||||||
|
20 moles = 0 # init target variable, (mol)
|
||||||
|
30 IF ((M <= 0) AND (sr_i < 1)) OR (sr_i = 1.0) THEN GOTO 310
|
||||||
|
40 sa = PARM(1) # reactive surface area, (m2)
|
||||||
|
|
||||||
|
100 r = 8.314462 # gas constant, (J K-1 mol-1)
|
||||||
|
110 dTi = (1 / TK) - (1 / 298.15) # (K-1)
|
||||||
|
120 ea_am = 30800 # activation energy am, (J mol-1)<---------change me
|
||||||
|
130 ea_nm = 30800 # activation energy nm, (J mol-1)<---------change me
|
||||||
|
rem ea_bm = 0 # activation energy bm, (J mol-1)<---------change me
|
||||||
|
150 log_k_am = -6.90 # reaction constant am<-----------------change me
|
||||||
|
160 log_k_nm = -7.90 # reaction constant nm<-----------------change me
|
||||||
|
rem log_k_bm = -99 # reaction constant bm<-------------------change me
|
||||||
|
180 n_am = 0.22 # H+ reaction order am<----------------------change me
|
||||||
|
rem n_bm = 0 # H+ reaction order bm<-----------------------change me
|
||||||
|
200 am = (10 ^ log_k_am) * EXP(-ea_am * dTi / r) * ACT("H+") ^ n_am
|
||||||
|
210 nm = (10 ^ log_k_nm) * EXP(-ea_nm * dTi / r)
|
||||||
|
rem bm = (10 ^ log_k_bm) * EXP(-ea_bm * dTi / r) * ACT("H+") ^ n_bm
|
||||||
|
|
||||||
|
300 moles = sa * (am + nm) * (1 - sr_i)
|
||||||
|
310 save moles * time
|
||||||
|
-end
|
||||||
|
|
||||||
|
END
|
||||||
80
poet/test/barite_het.pqi
Normal file
80
poet/test/barite_het.pqi
Normal file
@ -0,0 +1,80 @@
|
|||||||
|
## Initial: everywhere equilibrium with Celestite NB: The aqueous
|
||||||
|
## solution *resulting* from this calculation is to be used as initial
|
||||||
|
## state everywhere in the domain
|
||||||
|
SOLUTION 1
|
||||||
|
units mol/kgw
|
||||||
|
water 1
|
||||||
|
temperature 25
|
||||||
|
pH 7
|
||||||
|
pe 4
|
||||||
|
S(6) 1e-12
|
||||||
|
Sr 1e-12
|
||||||
|
Ba 1e-12
|
||||||
|
Cl 1e-12
|
||||||
|
PURE 1
|
||||||
|
Celestite 0.0 1
|
||||||
|
|
||||||
|
SAVE SOLUTION 2 # <- phreeqc keyword to store and later reuse these results
|
||||||
|
END
|
||||||
|
|
||||||
|
RUN_CELLS
|
||||||
|
-cells 1
|
||||||
|
|
||||||
|
COPY solution 1 2-3
|
||||||
|
|
||||||
|
## Here a 5x2 domain:
|
||||||
|
|
||||||
|
|---+---+---+---+---|
|
||||||
|
-> | 2 | 2 | 2 | 2 | 2 |
|
||||||
|
4 |---+---+---+---+---|
|
||||||
|
-> | 3 | 3 | 3 | 3 | 3 |
|
||||||
|
|---+---+---+---+---|
|
||||||
|
|
||||||
|
## East boundary: "injection" of solution 4. North, W, S boundaries: closed
|
||||||
|
|
||||||
|
## Here the two distinct zones: nr 2 with kinetics Celestite (initial
|
||||||
|
## amount is 0, is then allowed to precipitate) and nr 3 with kinetic
|
||||||
|
## Celestite and Barite (both initially > 0) where the actual
|
||||||
|
## replacement takes place
|
||||||
|
|
||||||
|
#USE SOLUTION 2 <- PHREEQC keyword to reuse the results from previous calculation
|
||||||
|
KINETICS 2
|
||||||
|
Celestite
|
||||||
|
-m 0 # Allowed to precipitate
|
||||||
|
-parms 10.0
|
||||||
|
-tol 1e-9
|
||||||
|
|
||||||
|
END
|
||||||
|
|
||||||
|
#USE SOLUTION 2 <- PHREEQC keyword to reuse the results from previous calculation
|
||||||
|
KINETICS 3
|
||||||
|
Barite
|
||||||
|
-m 0.001
|
||||||
|
-parms 50.
|
||||||
|
-tol 1e-9
|
||||||
|
Celestite
|
||||||
|
-m 1
|
||||||
|
-parms 10.0
|
||||||
|
-tol 1e-9
|
||||||
|
END
|
||||||
|
|
||||||
|
## A BaCl2 solution (nr 4) is "injected" from the left boundary:
|
||||||
|
SOLUTION 4
|
||||||
|
units mol/kgw
|
||||||
|
pH 7
|
||||||
|
water 1
|
||||||
|
temp 25
|
||||||
|
Ba 0.1
|
||||||
|
Cl 0.2
|
||||||
|
END
|
||||||
|
## NB: again, the *result* of the SOLUTION 4 script defines the
|
||||||
|
## concentration of all elements (+charge, tot H, tot O)
|
||||||
|
|
||||||
|
## Ideally, in the initial state SOLUTION 1 we should not have to
|
||||||
|
## define the 4 elemental concentrations (S(6), Sr, Ba and Cl) but
|
||||||
|
## obtain them having run once the scripts with the aqueous solution
|
||||||
|
## resulting from SOLUTION 1 once with KINETICS 2 and once with
|
||||||
|
## KINETICS 3.
|
||||||
|
|
||||||
|
RUN_CELLS
|
||||||
|
-cells 2-4
|
||||||
204
poet/test/testPhreeqcMatrix.cpp
Normal file
204
poet/test/testPhreeqcMatrix.cpp
Normal file
@ -0,0 +1,204 @@
|
|||||||
|
#include <cmath>
|
||||||
|
#include <gtest/gtest.h>
|
||||||
|
#include <vector>
|
||||||
|
|
||||||
|
#include <testPhreeqcMatrix.hpp>
|
||||||
|
|
||||||
|
#include "PhreeqcMatrix.hpp"
|
||||||
|
|
||||||
|
#define POET_TEST(name) TEST(TestPOET, name)
|
||||||
|
|
||||||
|
POET_TEST(PhreeqcInit) {
|
||||||
|
EXPECT_NO_THROW(
|
||||||
|
PhreeqcMatrix(base_test::phreeqc_database, base_test::script));
|
||||||
|
}
|
||||||
|
|
||||||
|
POET_TEST(PhreeqcMatrixOneSolution) {
|
||||||
|
PhreeqcMatrix pqc_mat(base_test::phreeqc_database, base_test::script);
|
||||||
|
const auto ids = pqc_mat.getIds();
|
||||||
|
EXPECT_EQ(ids.size(), 1);
|
||||||
|
EXPECT_EQ(ids[0], 1);
|
||||||
|
|
||||||
|
PhreeqcMatrix::STLExport exported_init = pqc_mat.get();
|
||||||
|
// ID + H,O,Charge + 4 Solutions + 4 Equil incl. params
|
||||||
|
EXPECT_EQ(exported_init.names.size(), 12);
|
||||||
|
|
||||||
|
EXPECT_EQ(exported_init.names, base_test::expected_names);
|
||||||
|
for (std::size_t i = 0; i < exported_init.values.size(); ++i) {
|
||||||
|
EXPECT_NEAR(exported_init.values[i], base_test::expected_values[i],
|
||||||
|
base_test::expected_errors[i]);
|
||||||
|
}
|
||||||
|
|
||||||
|
auto dumps = pqc_mat.getDumpStringsPQI();
|
||||||
|
EXPECT_EQ(dumps.size(), 1);
|
||||||
|
EXPECT_GT(dumps[1].size(), 1);
|
||||||
|
|
||||||
|
const auto kinetics = pqc_mat.getKineticsNames(1);
|
||||||
|
EXPECT_EQ(kinetics.size(), 0);
|
||||||
|
|
||||||
|
const auto equilibrium = pqc_mat.getEquilibriumNames(1);
|
||||||
|
EXPECT_EQ(equilibrium.size(), 2);
|
||||||
|
|
||||||
|
const auto expected_equilibrium =
|
||||||
|
std::vector<std::string>({"Calcite", "Dolomite"});
|
||||||
|
EXPECT_EQ(equilibrium, expected_equilibrium);
|
||||||
|
}
|
||||||
|
|
||||||
|
POET_TEST(PhreeqcMatrixMultiSolution) {
|
||||||
|
PhreeqcMatrix pqc_mat(barite_test::database, barite_test::script);
|
||||||
|
|
||||||
|
const auto ids = pqc_mat.getIds();
|
||||||
|
EXPECT_EQ(ids.size(), 4);
|
||||||
|
EXPECT_EQ(ids[0], 1);
|
||||||
|
EXPECT_EQ(ids[1], 2);
|
||||||
|
EXPECT_EQ(ids[2], 3);
|
||||||
|
EXPECT_EQ(ids[3], 4);
|
||||||
|
|
||||||
|
PhreeqcMatrix::STLExport exported = pqc_mat.get();
|
||||||
|
|
||||||
|
EXPECT_EQ(exported.names, barite_test::expected_names);
|
||||||
|
for (std::size_t i = 0; i < exported.names.size(); i++) {
|
||||||
|
if (i > 8 && i < 13) {
|
||||||
|
EXPECT_TRUE(std::isnan(exported.values[i]));
|
||||||
|
continue;
|
||||||
|
}
|
||||||
|
EXPECT_NEAR(exported.values[i], barite_test::expected_values_line_one[i],
|
||||||
|
barite_test::expected_errors[i]);
|
||||||
|
}
|
||||||
|
|
||||||
|
auto dumps = pqc_mat.getDumpStringsPQI();
|
||||||
|
EXPECT_EQ(dumps.size(), 4);
|
||||||
|
|
||||||
|
const auto kinetics_sol_1 = pqc_mat.getKineticsNames(1);
|
||||||
|
EXPECT_EQ(kinetics_sol_1.size(), 0);
|
||||||
|
|
||||||
|
const auto equilibrium_sol_1 = pqc_mat.getEquilibriumNames(1);
|
||||||
|
const auto expected_equilibrium_sol_1 =
|
||||||
|
std::vector<std::string>({"Celestite"});
|
||||||
|
EXPECT_EQ(equilibrium_sol_1, expected_equilibrium_sol_1);
|
||||||
|
|
||||||
|
const auto equilibrium_sol_2 = pqc_mat.getEquilibriumNames(2);
|
||||||
|
EXPECT_EQ(equilibrium_sol_2.size(), 0);
|
||||||
|
|
||||||
|
const auto kinetics_sol_2 = pqc_mat.getKineticsNames(2);
|
||||||
|
const auto expected_kinetics_sol_2 = std::vector<std::string>({"Celestite"});
|
||||||
|
EXPECT_EQ(kinetics_sol_2, expected_kinetics_sol_2);
|
||||||
|
|
||||||
|
const auto kinetics_sol_3 = pqc_mat.getKineticsNames(3);
|
||||||
|
const auto expected_kinetics_sol_3 =
|
||||||
|
std::vector<std::string>({"Barite", "Celestite"});
|
||||||
|
EXPECT_EQ(kinetics_sol_3, expected_kinetics_sol_3);
|
||||||
|
|
||||||
|
EXPECT_EQ(pqc_mat.getKineticsNames(4).size(), 0);
|
||||||
|
EXPECT_EQ(pqc_mat.getEquilibriumNames(4).size(), 0);
|
||||||
|
}
|
||||||
|
|
||||||
|
POET_TEST(PhreeqcMatrixCtor) {
|
||||||
|
PhreeqcMatrix pqc_mat(barite_test::database, barite_test::script);
|
||||||
|
PhreeqcMatrix pqc_mat_copy(pqc_mat);
|
||||||
|
PhreeqcMatrix pqc_mat_move(std::move(pqc_mat_copy));
|
||||||
|
|
||||||
|
const auto ids = pqc_mat_move.getIds();
|
||||||
|
EXPECT_EQ(ids.size(), 4);
|
||||||
|
EXPECT_EQ(ids[0], 1);
|
||||||
|
EXPECT_EQ(ids[1], 2);
|
||||||
|
EXPECT_EQ(ids[2], 3);
|
||||||
|
EXPECT_EQ(ids[3], 4);
|
||||||
|
|
||||||
|
PhreeqcMatrix::STLExport exported = pqc_mat_move.get();
|
||||||
|
|
||||||
|
EXPECT_EQ(exported.names, barite_test::expected_names);
|
||||||
|
for (std::size_t i = 0; i < exported.names.size(); i++) {
|
||||||
|
if (i > 8 && i < 13) {
|
||||||
|
EXPECT_TRUE(std::isnan(exported.values[i]));
|
||||||
|
continue;
|
||||||
|
}
|
||||||
|
EXPECT_NEAR(exported.values[i], barite_test::expected_values_line_one[i],
|
||||||
|
barite_test::expected_errors[i]);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
POET_TEST(PhreeqcMatrixOperator) {
|
||||||
|
PhreeqcMatrix pqc_mat(barite_test::database, barite_test::script);
|
||||||
|
PhreeqcMatrix pqc_mat_copy = pqc_mat;
|
||||||
|
|
||||||
|
const auto ids = pqc_mat_copy.getIds();
|
||||||
|
|
||||||
|
EXPECT_EQ(ids.size(), 4);
|
||||||
|
EXPECT_EQ(ids[0], 1);
|
||||||
|
EXPECT_EQ(ids[1], 2);
|
||||||
|
EXPECT_EQ(ids[2], 3);
|
||||||
|
EXPECT_EQ(ids[3], 4);
|
||||||
|
|
||||||
|
PhreeqcMatrix::STLExport exported = pqc_mat_copy.get();
|
||||||
|
|
||||||
|
EXPECT_EQ(exported.names, barite_test::expected_names);
|
||||||
|
for (std::size_t i = 0; i < exported.names.size(); i++) {
|
||||||
|
if (i > 8 && i < 13) {
|
||||||
|
EXPECT_TRUE(std::isnan(exported.values[i]));
|
||||||
|
continue;
|
||||||
|
}
|
||||||
|
EXPECT_NEAR(exported.values[i], barite_test::expected_values_line_one[i],
|
||||||
|
barite_test::expected_errors[i]);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
POET_TEST(PhreeqcMatrixRvalueManipulation) {
|
||||||
|
PhreeqcMatrix pqc_mat(barite_test::database, barite_test::script);
|
||||||
|
|
||||||
|
PhreeqcMatrix pqc_erased = pqc_mat.erase({1});
|
||||||
|
|
||||||
|
const std::vector<int> expected_ids_erased = {2, 3, 4};
|
||||||
|
|
||||||
|
PhreeqcMatrix::STLExport exported_erased = pqc_erased.get();
|
||||||
|
|
||||||
|
EXPECT_EQ(pqc_erased.getIds(), expected_ids_erased);
|
||||||
|
EXPECT_EQ(exported_erased.names, barite_test::expected_names_erased);
|
||||||
|
EXPECT_EQ(exported_erased.values.size(),
|
||||||
|
exported_erased.names.size() * expected_ids_erased.size());
|
||||||
|
|
||||||
|
PhreeqcMatrix pqc_mat_subset = pqc_mat.subset({1});
|
||||||
|
|
||||||
|
const std::vector<int> expected_ids_subset = {1};
|
||||||
|
|
||||||
|
PhreeqcMatrix::STLExport exported_subset = pqc_mat_subset.get();
|
||||||
|
|
||||||
|
EXPECT_EQ(pqc_mat_subset.getIds(), expected_ids_subset);
|
||||||
|
EXPECT_EQ(exported_subset.names, barite_test::expected_names_subset);
|
||||||
|
|
||||||
|
pqc_mat = pqc_mat.subset({1});
|
||||||
|
|
||||||
|
exported_subset = pqc_mat.get();
|
||||||
|
|
||||||
|
EXPECT_EQ(pqc_mat_subset.getIds(), expected_ids_subset);
|
||||||
|
EXPECT_EQ(exported_subset.names, barite_test::expected_names_subset);
|
||||||
|
}
|
||||||
|
|
||||||
|
POET_TEST(PhreeqcMatrixColumnMajorExport) {
|
||||||
|
PhreeqcMatrix pqc_mat(barite_test::database, barite_test::script);
|
||||||
|
|
||||||
|
pqc_mat = pqc_mat.subset({2, 3});
|
||||||
|
|
||||||
|
PhreeqcMatrix::STLExport exported =
|
||||||
|
pqc_mat.get(PhreeqcMatrix::VectorExportType::COLUMN_MAJOR);
|
||||||
|
|
||||||
|
const auto ids = pqc_mat.getIds();
|
||||||
|
EXPECT_EQ(ids.size(), 2);
|
||||||
|
EXPECT_EQ(ids[0], 2);
|
||||||
|
EXPECT_EQ(ids[1], 3);
|
||||||
|
|
||||||
|
EXPECT_EQ(exported.names, barite_test::expected_names_erased);
|
||||||
|
EXPECT_EQ(exported.values.size(), exported.names.size() * ids.size());
|
||||||
|
|
||||||
|
EXPECT_EQ(exported.values[0], 2);
|
||||||
|
EXPECT_EQ(exported.values[1], 3);
|
||||||
|
}
|
||||||
|
|
||||||
|
POET_TEST(PhreeqcMatrixBracketOperator) {
|
||||||
|
PhreeqcMatrix pqc_mat(base_test::phreeqc_database, base_test::script);
|
||||||
|
|
||||||
|
EXPECT_NO_THROW(pqc_mat(1, "H"));
|
||||||
|
EXPECT_NEAR(pqc_mat(1, "H"), base_test::expected_values[1], 1e-5);
|
||||||
|
EXPECT_ANY_THROW(pqc_mat(1, "J"));
|
||||||
|
EXPECT_ANY_THROW(pqc_mat(2, "H"));
|
||||||
|
}
|
||||||
105
poet/test/testPhreeqcMatrix.hpp.in
Normal file
105
poet/test/testPhreeqcMatrix.hpp.in
Normal file
@ -0,0 +1,105 @@
|
|||||||
|
#pragma once
|
||||||
|
|
||||||
|
#include <limits>
|
||||||
|
#include <string>
|
||||||
|
#include <vector>
|
||||||
|
|
||||||
|
namespace base_test {
|
||||||
|
const std::string script = R"(SOLUTION 1
|
||||||
|
units mol/kgw
|
||||||
|
temp 25
|
||||||
|
Ca 0.1
|
||||||
|
Mg 0.1
|
||||||
|
Cl 0.5 charge
|
||||||
|
Na 0.1
|
||||||
|
PURE 1
|
||||||
|
Calcite 0.0 1
|
||||||
|
Dolomite 0.0 0
|
||||||
|
## RUN_CELLS
|
||||||
|
## -cells 1
|
||||||
|
END)";
|
||||||
|
|
||||||
|
const std::vector<std::string> expected_names = {
|
||||||
|
"ID", "H", "O", "Charge", "Ca", "Cl",
|
||||||
|
"Mg", "Na", "Calcite_eq", "Calcite_si", "Dolomite_eq", "Dolomite_si"};
|
||||||
|
|
||||||
|
const std::vector<double> expected_values = {1,
|
||||||
|
111.01243521533338,
|
||||||
|
55.506218386370165,
|
||||||
|
-4.7941959748097226e-13,
|
||||||
|
0.1,
|
||||||
|
0.5,
|
||||||
|
0.1,
|
||||||
|
0.1,
|
||||||
|
1,
|
||||||
|
0,
|
||||||
|
0,
|
||||||
|
0};
|
||||||
|
const std::vector<double> expected_errors = {
|
||||||
|
0, 1e-3, 1e-3, 1e-15, 1e-5, 1e-5, 1e-5, 1e-5, 1e-5, 0, 1e-5, 0};
|
||||||
|
|
||||||
|
const std::string phreeqc_database = R"database(@POET_PHREEQCDAT_DB@)database";
|
||||||
|
} // namespace base_test
|
||||||
|
|
||||||
|
namespace barite_test {
|
||||||
|
const std::string script = R"barite(@POET_BARITE_PQI@)barite";
|
||||||
|
const std::string database = R"barite(@POET_BARITE_DB@)barite";
|
||||||
|
|
||||||
|
const std::vector<std::string> expected_names = {"ID",
|
||||||
|
"H",
|
||||||
|
"O",
|
||||||
|
"Charge",
|
||||||
|
"Ba",
|
||||||
|
"Cl",
|
||||||
|
"S(-2)",
|
||||||
|
"S(6)",
|
||||||
|
"Sr",
|
||||||
|
"Barite",
|
||||||
|
"Barite_p1",
|
||||||
|
"Celestite",
|
||||||
|
"Celestite_p1",
|
||||||
|
"Celestite_eq",
|
||||||
|
"Celestite_si"};
|
||||||
|
|
||||||
|
const std::vector<double> expected_values_line_one = {
|
||||||
|
1,
|
||||||
|
111.01243359383071,
|
||||||
|
55.508698688362124,
|
||||||
|
-1.2153078399577636e-09,
|
||||||
|
1.0000001848805677e-12,
|
||||||
|
1.0000000116187218e-12,
|
||||||
|
0,
|
||||||
|
0.00062047270964839664,
|
||||||
|
0.00062047270964840271,
|
||||||
|
std::numeric_limits<double>::quiet_NaN(),
|
||||||
|
std::numeric_limits<double>::quiet_NaN(),
|
||||||
|
std::numeric_limits<double>::quiet_NaN(),
|
||||||
|
std::numeric_limits<double>::quiet_NaN(),
|
||||||
|
0.99937952729135193,
|
||||||
|
0};
|
||||||
|
|
||||||
|
const std::vector<double> expected_errors = {
|
||||||
|
0,
|
||||||
|
1e-3,
|
||||||
|
1e-3,
|
||||||
|
1e-15,
|
||||||
|
1e-5,
|
||||||
|
1e-5,
|
||||||
|
1e-5,
|
||||||
|
1e-5,
|
||||||
|
1e-5,
|
||||||
|
std::numeric_limits<double>::quiet_NaN(),
|
||||||
|
std::numeric_limits<double>::quiet_NaN(),
|
||||||
|
std::numeric_limits<double>::quiet_NaN(),
|
||||||
|
std::numeric_limits<double>::quiet_NaN(),
|
||||||
|
1e-5,
|
||||||
|
0};
|
||||||
|
|
||||||
|
const std::vector<std::string> expected_names_erased = {
|
||||||
|
"ID", "H", "O", "Charge", "Ba", "Cl", "S(-2)",
|
||||||
|
"S(6)", "Sr", "Barite", "Barite_p1", "Celestite", "Celestite_p1"};
|
||||||
|
|
||||||
|
const std::vector<std::string> expected_names_subset = {
|
||||||
|
"ID", "H", "O", "Charge", "Ba", "Cl",
|
||||||
|
"S(-2)", "S(6)", "Sr", "Celestite_eq", "Celestite_si"};
|
||||||
|
} // namespace barite_test
|
||||||
@ -888,6 +888,8 @@ public:
|
|||||||
virtual bool output_open(const char *file_name, std::ios_base::openmode mode = std::ios_base::out);
|
virtual bool output_open(const char *file_name, std::ios_base::openmode mode = std::ios_base::out);
|
||||||
virtual bool punch_open(const char *file_name, std::ios_base::openmode mode = std::ios_base::out, int n_user = 1);
|
virtual bool punch_open(const char *file_name, std::ios_base::openmode mode = std::ios_base::out, int n_user = 1);
|
||||||
|
|
||||||
|
Phreeqc *GetPhreeqcPtr(void) { return this->PhreeqcPtr; };
|
||||||
|
|
||||||
protected:
|
protected:
|
||||||
int EndRow(void);
|
int EndRow(void);
|
||||||
void AddSelectedOutput(const char* name, const char* format, va_list argptr);
|
void AddSelectedOutput(const char* name, const char* format, va_list argptr);
|
||||||
|
|||||||
75
src/phreeqcpp/GFZ.cpp
Normal file
75
src/phreeqcpp/GFZ.cpp
Normal file
@ -0,0 +1,75 @@
|
|||||||
|
#include "Phreeqc.h"
|
||||||
|
|
||||||
|
std::vector<std::string> Phreeqc::find_all_valence_states(
|
||||||
|
const std::vector<std::string> &&solution_names, const std::size_t offset) {
|
||||||
|
std::vector<std::string> solution_with_valences(
|
||||||
|
solution_names.begin(), solution_names.begin() + offset);
|
||||||
|
|
||||||
|
// to avoid duplicates store already evaluated master species
|
||||||
|
std::set<std::string> master_species_found;
|
||||||
|
|
||||||
|
for (std::size_t i = offset; i < solution_names.size(); i++) {
|
||||||
|
const auto master_primary =
|
||||||
|
master_bsearch_primary(solution_names[i].c_str());
|
||||||
|
|
||||||
|
assert(master_primary != NULL);
|
||||||
|
|
||||||
|
const std::string master_species(master_primary->elt->name);
|
||||||
|
|
||||||
|
// in case we already found valences for this master species we skip it
|
||||||
|
if (master_species_found.find(master_species) !=
|
||||||
|
master_species_found.end()) {
|
||||||
|
continue;
|
||||||
|
}
|
||||||
|
|
||||||
|
master_species_found.insert(master_species);
|
||||||
|
|
||||||
|
bool has_valences = false;
|
||||||
|
std::size_t inner_loop_j = 0;
|
||||||
|
std::size_t last_valence = inner_loop_j;
|
||||||
|
|
||||||
|
// we HAVE to assume master species are already sorted!
|
||||||
|
// loop over all known master species
|
||||||
|
for (inner_loop_j = 0; inner_loop_j < master.size(); inner_loop_j++) {
|
||||||
|
std::string curr_master_species(master[inner_loop_j]->elt->name);
|
||||||
|
|
||||||
|
if (curr_master_species == master_species) {
|
||||||
|
last_valence = inner_loop_j;
|
||||||
|
|
||||||
|
while (last_valence < master.size() - 1) {
|
||||||
|
std::string next_master_species(master[last_valence + 1]->elt->name);
|
||||||
|
|
||||||
|
// check if the next name starts with the current master species
|
||||||
|
if (next_master_species.compare(0, master_species.size(),
|
||||||
|
master_species) != 0) {
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
|
||||||
|
// check if the next character is an opening parenthesis
|
||||||
|
if (next_master_species[master_species.size()] != '(') {
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
|
||||||
|
// if this is all true we have a valence state
|
||||||
|
has_valences = true;
|
||||||
|
last_valence++;
|
||||||
|
}
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// in case we found valences we add them to the solution
|
||||||
|
if (has_valences) {
|
||||||
|
// skip the master species and only add the valences
|
||||||
|
for (std::size_t k = inner_loop_j + 1; k <= last_valence; k++) {
|
||||||
|
solution_with_valences.push_back(master[k]->elt->name);
|
||||||
|
}
|
||||||
|
continue;
|
||||||
|
}
|
||||||
|
|
||||||
|
// otherwise we just add the master species without any valences
|
||||||
|
solution_with_valences.push_back(master_species);
|
||||||
|
}
|
||||||
|
|
||||||
|
return solution_with_valences;
|
||||||
|
}
|
||||||
@ -1833,14 +1833,20 @@ protected:
|
|||||||
friend class IPhreeqcMMS;
|
friend class IPhreeqcMMS;
|
||||||
friend class IPhreeqcPhast;
|
friend class IPhreeqcPhast;
|
||||||
friend class PhreeqcRM;
|
friend class PhreeqcRM;
|
||||||
friend class PhreeqcInit;
|
friend class PhreeqcInit;
|
||||||
friend class PhreeqcEngine;
|
friend class PhreeqcEngine;
|
||||||
|
|
||||||
std::vector<int> keycount; // used to mark keywords that have been read
|
std::vector<int> keycount; // used to mark keywords that have been read
|
||||||
|
|
||||||
|
int lol;
|
||||||
|
|
||||||
public:
|
public:
|
||||||
static const class const_iso iso_defaults[];
|
std::vector<std::string>
|
||||||
static const int count_iso_defaults;
|
find_all_valence_states(const std::vector<std::string> &&solution_names,
|
||||||
|
const std::size_t offset);
|
||||||
|
|
||||||
|
static const class const_iso iso_defaults[];
|
||||||
|
static const int count_iso_defaults;
|
||||||
};
|
};
|
||||||
#endif /* _INC_PHREEQC_H */
|
#endif /* _INC_PHREEQC_H */
|
||||||
|
|
||||||
|
|||||||
Loading…
x
Reference in New Issue
Block a user