mirror of
https://git.gfz-potsdam.de/naaice/iphreeqc.git
synced 2025-12-15 16:18:22 +01:00
Improve usability outside POET
This commit is contained in:
Max Luebke
Max Lübke
parent
e02969c891
commit
ee940d20e4
148
.gitlab-ci.yml
148
.gitlab-ci.yml
@ -1,134 +1,22 @@
|
||||
#
|
||||
# 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
|
||||
image: ubuntu:24.04
|
||||
|
||||
before_script:
|
||||
- eval $(ssh-agent -s)
|
||||
- echo "${SSH_PRIVATE_KEY_ENC}" | base64 --decode | tr -d '\r' | ssh-add -
|
||||
- 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:
|
||||
stage: sync
|
||||
- apt-get update -y
|
||||
- apt-get install -y build-essential cmake git
|
||||
|
||||
##
|
||||
## 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
|
||||
|
||||
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}
|
||||
stages:
|
||||
- test
|
||||
|
||||
##
|
||||
## 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
|
||||
test:
|
||||
stage: test
|
||||
script:
|
||||
- mkdir _build && cd _build
|
||||
- cmake -DBUILD_SHARED_LIBS=ON -DCMAKE_BUILD_TYPE=Release ..
|
||||
- make -j$(nproc)
|
||||
- ctest --output-junit test_results.xml
|
||||
artifacts:
|
||||
when: always
|
||||
paths:
|
||||
- _build/test_results.xml
|
||||
reports:
|
||||
junit: _build/test_results.xml
|
||||
@ -220,6 +220,7 @@ target_sources(IPhreeqc
|
||||
src/phreeqcpp/UserPunch.cpp
|
||||
src/phreeqcpp/UserPunch.h
|
||||
src/phreeqcpp/utilities.cpp
|
||||
src/phreeqcpp/GFZ.cpp
|
||||
src/thread.h
|
||||
src/Var.c
|
||||
src/Var.h
|
||||
@ -392,9 +393,6 @@ if (BUILD_CLR_LIBS)
|
||||
|
||||
endif()
|
||||
|
||||
add_subdirectory(poet)
|
||||
add_subdirectory(app)
|
||||
|
||||
include (CTest)
|
||||
|
||||
if (STANDALONE_BUILD)
|
||||
@ -459,6 +457,9 @@ if (STANDALONE_BUILD)
|
||||
endif()
|
||||
endif()
|
||||
|
||||
add_subdirectory(poet)
|
||||
# add_subdirectory(app)
|
||||
|
||||
# get_cmake_property(_variableNames VARIABLES)
|
||||
# list (SORT _variableNames)
|
||||
# foreach (_variableName ${_variableNames})
|
||||
|
||||
@ -5,8 +5,11 @@ add_library(IPhreeqcPOET ${PQC_POET_SRC})
|
||||
target_link_libraries(IPhreeqcPOET PUBLIC IPhreeqc)
|
||||
target_include_directories(IPhreeqcPOET PUBLIC
|
||||
$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
|
||||
$<INSTALL_INTERFACE:include>
|
||||
PRIVATE src)
|
||||
$<INSTALL_INTERFACE:include>)
|
||||
|
||||
# Set C++23 standard
|
||||
target_compile_features(IPhreeqcPOET PUBLIC cxx_std_20)
|
||||
# Set C++20 standard
|
||||
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
|
||||
|
||||
#include "POETInit.hpp"
|
||||
|
||||
#include "PhreeqcMatrix.hpp"
|
||||
#include <memory>
|
||||
#include <vector>
|
||||
|
||||
/**
|
||||
* @brief Class for running Phreeqc wrappped in POET
|
||||
*
|
||||
* Direct interface to Phreeqc, without utilizing *_MODIFY keywords/scripts
|
||||
to
|
||||
* set new values. Use with already initialized Phreeqc config.
|
||||
* Direct interface to Phreeqc, without utilizing Phreeqc's internal parser
|
||||
* to set/get values. To start a simulation, the interpreter is used internally.
|
||||
* No need for the user to directly interact with the interpreter.
|
||||
*
|
||||
*/
|
||||
class PhreeqcEngine {
|
||||
@ -19,11 +18,13 @@ public:
|
||||
* @brief Construct a new Phreeqc Engine object
|
||||
*
|
||||
* 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
|
||||
@ -35,8 +36,8 @@ public:
|
||||
* @brief Siimulate a cell for a given time step
|
||||
*
|
||||
* @param cell_values Vector containing the input values for the cell
|
||||
* (**including the ID**). Output values are written back in place to this
|
||||
* vector!
|
||||
* (*including the ID*). Output values are written back in place to this
|
||||
* vector.
|
||||
* @param time_step Time step to simulate in seconds
|
||||
*/
|
||||
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 <cstddef>
|
||||
#include <iomanip>
|
||||
#include <regex>
|
||||
#include <span>
|
||||
#include <sstream>
|
||||
|
||||
#include <IPhreeqc.hpp>
|
||||
#include <Phreeqc.h>
|
||||
#include <string>
|
||||
#include <vector>
|
||||
|
||||
#include "PhreeqcMatrix.hpp"
|
||||
#include "Wrapper/EquilibriumWrapper.hpp"
|
||||
#include "Wrapper/ExchangeWrapper.hpp"
|
||||
#include "Wrapper/KineticWrapper.hpp"
|
||||
@ -14,8 +18,7 @@
|
||||
#include "Wrapper/SurfaceWrapper.hpp"
|
||||
class PhreeqcEngine::Impl : public IPhreeqc {
|
||||
public:
|
||||
Impl(const POETConfig &config);
|
||||
void init_wrappers(const POETInitCell &cell);
|
||||
Impl(const PhreeqcMatrix &pqc_mat, const int cell_id);
|
||||
void run(double time_step);
|
||||
|
||||
cxxSolution *Get_solution(std::size_t n) {
|
||||
@ -53,20 +56,48 @@ public:
|
||||
bool has_kinetics = false;
|
||||
bool has_equilibrium = 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)
|
||||
: impl(std::make_unique<Impl>(config)) {}
|
||||
|
||||
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(const PhreeqcMatrix &pqc_mat, const int cell_id)
|
||||
: impl(std::make_unique<Impl>(pqc_mat, cell_id)) {}
|
||||
|
||||
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,
|
||||
double time_step) {
|
||||
// skip ID
|
||||
@ -86,7 +117,7 @@ void PhreeqcEngine::Impl::run(double time_step) {
|
||||
this->RunString(runs_string.c_str());
|
||||
}
|
||||
|
||||
void PhreeqcEngine::Impl::init_wrappers(const POETInitCell &cell) {
|
||||
void PhreeqcEngine::Impl::init_wrappers(const InitCell &cell) {
|
||||
|
||||
// Solutions
|
||||
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;
|
||||
|
||||
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");
|
||||
|
||||
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 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:
|
||||
int EndRow(void);
|
||||
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 IPhreeqcPhast;
|
||||
friend class PhreeqcRM;
|
||||
friend class PhreeqcInit;
|
||||
friend class PhreeqcEngine;
|
||||
friend class PhreeqcInit;
|
||||
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:
|
||||
static const class const_iso iso_defaults[];
|
||||
static const int count_iso_defaults;
|
||||
std::vector<std::string>
|
||||
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 */
|
||||
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user