diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 6267b622..a64ef30e 100644 --- a/.gitlab-ci.yml +++ b/.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 \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index f3a42f2a..8ffd67fd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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 @@ -398,9 +399,6 @@ if (BUILD_CLR_LIBS) endif() -add_subdirectory(poet) -add_subdirectory(app) - include (CTest) if (STANDALONE_BUILD) @@ -472,6 +470,9 @@ if (STANDALONE_BUILD) endif() endif() +add_subdirectory(poet) +# add_subdirectory(app) + # get_cmake_property(_variableNames VARIABLES) # list (SORT _variableNames) # foreach (_variableName ${_variableNames}) diff --git a/poet/CMakeLists.txt b/poet/CMakeLists.txt index f6f4b0d6..4d734adb 100644 --- a/poet/CMakeLists.txt +++ b/poet/CMakeLists.txt @@ -5,8 +5,11 @@ add_library(IPhreeqcPOET ${PQC_POET_SRC}) target_link_libraries(IPhreeqcPOET PUBLIC IPhreeqc) target_include_directories(IPhreeqcPOET PUBLIC $ - $ - PRIVATE src) + $) -# Set C++23 standard -target_compile_features(IPhreeqcPOET PUBLIC cxx_std_20) \ No newline at end of file +# Set C++20 standard +target_compile_features(IPhreeqcPOET PUBLIC cxx_std_20) + +if (BUILD_TESTING AND STANDALONE_BUILD) + add_subdirectory(test) +endif() \ No newline at end of file diff --git a/poet/include/POETInit.hpp b/poet/include/POETInit.hpp deleted file mode 100644 index 2bdc12e1..00000000 --- a/poet/include/POETInit.hpp +++ /dev/null @@ -1,20 +0,0 @@ -#pragma once - -#include -#include - -struct POETInitCell { - std::vector solutions; - std::vector solution_primaries; - std::vector exchanger; - std::vector kinetics; - std::vector equilibrium; - std::vector surface_comps; - std::vector surface_charges; -}; - -struct POETConfig { - std::string database; - std::string input_script; - POETInitCell cell; -}; \ No newline at end of file diff --git a/poet/include/PhreeqcEngine.hpp b/poet/include/PhreeqcEngine.hpp index 91ae9a2c..5bb3132c 100644 --- a/poet/include/PhreeqcEngine.hpp +++ b/poet/include/PhreeqcEngine.hpp @@ -1,16 +1,15 @@ #pragma once -#include "POETInit.hpp" - +#include "PhreeqcMatrix.hpp" #include #include /** * @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 &cell_values, double time_step); diff --git a/poet/include/PhreeqcInit.hpp b/poet/include/PhreeqcInit.hpp deleted file mode 100644 index 02433678..00000000 --- a/poet/include/PhreeqcInit.hpp +++ /dev/null @@ -1,247 +0,0 @@ -#pragma once - -#include - -#include -#include -#include -#include -#include -#include -#include -#include - -#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 names; - std::vector ids; - std::vector> 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 Contains the raw input scripts for each - * cell - */ - std::map raw_dumps(); - - using essential_names = std::array, 5>; - using ModulesArray = std::array; - - /** - * @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 &cell_ids); - - /** - * @brief Get solution primaries found by PhreeqcInit - * - * @return std::vector Vector containing only the solution - * primary names of the solution defined by the input script - */ - std::vector getSolutionPrimaries() { - return std::vector(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 Whole vector of solution names for the - * cell - */ - std::vector 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 Whole vector of exchange names for the - * cell. Empty if no exchange is defined. - */ - std::vector 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 Whole vector of kinetics names for the - * cell. Empty if no kinetics are defined. - */ - std::vector 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 Whole vector of equilibrium names for the - * cell. Empty if no equilibrium is defined. - */ - std::vector 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 Whole vector of surface component names - * for the cell. Empty if no surface is defined. - */ - std::vector 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 Whole vector of surface charge names for - * the cell. Empty if no surface is defined. - */ - std::vector 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 dump_solution_names(int cell_number); - void dump_reactant_names(int cell_number, - const std::vector union_sol_names, - essential_names &names); - - std::vector - find_all_valence_states(const std::vector &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> exchanger; - std::map> kinetics; - std::map> equilibrium; - std::map> surface_comps; - std::map> surface_charge; - - bool contains_id(const std::map> &map, - int id) const { -#if __cplusplus >= 202002L - return map.contains(id); -#else - return map.find(id) != map.end(); -#endif - } - - std::set surface_primaries; - - using RawMap = std::map>>; - - essential_names union_raws(const RawMap &raws); - - std::vector> - 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 &values); - - std::string subExchangeName(std::string name); - - std::vector - get_essential_values_init(std::size_t cell_number, - const std::vector &order); - - RawMap raw_initials; -}; \ No newline at end of file diff --git a/poet/include/PhreeqcMatrix.hpp b/poet/include/PhreeqcMatrix.hpp new file mode 100644 index 00000000..3cf8eaf6 --- /dev/null +++ b/poet/include/PhreeqcMatrix.hpp @@ -0,0 +1,274 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +#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 &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 &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 names; + std::vector 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(PhreeqcComponent::EXCHANGE), + KINETICS, + EQUILIBRIUM, + SURACE_COMP, + SURFACE_CHARGE + } type; + + std::string name; + }; + + /** + * @brief Get all found solution names + * + * @return std::vector Vector containing all solution names. + */ + std::vector + getSolutionNames(bool include_h_o_charge = false) const; + + /** + * @brief Get solution total names of all found solutions (excluding H, O, + * Charge) + * + * @return std::vector Names of all found solutions (excluding H, + * O, Charge) + */ + std::vector 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 Whole vector of exchange names for the + * cell. Empty if no exchange is defined. + */ + std::vector 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 Whole vector of kinetics names for the + * cell. Empty if no kinetics are defined. + */ + std::vector 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 Whole vector of equilibrium names for the + * cell. Empty if no equilibrium is defined. + */ + std::vector 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 Whole vector of surface component names + * for the cell. Empty if no surface is defined. + */ + std::vector 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 Whole vector of surface charge names for + * the cell. Empty if no surface is defined. + */ + std::vector getSurfaceChargeNames(int cell_id) const; + + /** + * @brief Get all cell IDs stored in the PhreeqcMatrix. + * + * @return std::vector IDs of all cells stored in the PhreeqcMatrix. + */ + std::vector getIds() const; + + // std::array getComponentCount(int cell_id) const; + + /** + * @brief Dump all cells into a **DUMP** format of Phreeqc. + * + * @return std::map Map containing the cell ID as key and + * the exported DUMP string as value. + */ + std::map 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> _m_map; + std::map> _m_internal_names; + + std::set _m_surface_primaries; + + void initialize(); + + void remove_NaNs(); + + std::shared_ptr _m_pqc; + std::string _m_database; +}; \ No newline at end of file diff --git a/poet/src/Engine.cpp b/poet/src/Engine.cpp index c74b73da..8fad4002 100644 --- a/poet/src/Engine.cpp +++ b/poet/src/Engine.cpp @@ -1,12 +1,16 @@ #include "PhreeqcEngine.hpp" #include #include +#include #include #include #include #include +#include +#include +#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 solutions; + std::vector exchanger; + std::vector kinetics; + std::vector equilibrium; + std::vector surface_comps; + std::vector surface_charges; + std::vector solution_primaries; + }; + void init_wrappers(const InitCell &cell); }; -PhreeqcEngine::PhreeqcEngine(const POETConfig &config) - : impl(std::make_unique(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(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 &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 = diff --git a/poet/src/Extraction.cpp b/poet/src/Extraction.cpp deleted file mode 100644 index d8b1ea1f..00000000 --- a/poet/src/Extraction.cpp +++ /dev/null @@ -1,204 +0,0 @@ -#include - -#include "Solution.h" - -#include -#include -#include -#include - -static inline std::vector -unionStringVectors(const std::vector &a, - const std::vector &b) { - std::vector result; - - std::set_union(a.begin(), a.end(), b.begin(), b.end(), - std::back_inserter(result)); - - return result; -} - -static std::vector -createConcVector(const std::vector &conc_names, - const std::vector::const_iterator &conc_it, - const std::vector &new_names) { - std::vector 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::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 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 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 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> -PhreeqcInit::conc_from_essentials(const RawMap &raws, - const essential_names &names) { - std::vector> values; - - for (const auto &[key, val] : raws) { - std::vector combined_conc_vec; - - for (std::size_t i = 0, offset = 0; i < names.size(); i++) { - std::vector 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> values; - - const PhreeqcInit::essential_names ess_names = - this->union_raws(this->raw_initials); - - const std::vector> conc_values = - this->conc_from_essentials(this->raw_initials, ess_names); - - std::vector conc_names; - - for (const auto &i : ess_names) { - conc_names.insert(conc_names.end(), i.begin(), i.end()); - } - - std::vector 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 &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 &use_ids) { -// std::vector> 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> conc_values = -// this->conc_from_essentials(raws, ess_names); - -// std::vector 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}; -// } \ No newline at end of file diff --git a/poet/src/InitGetSet.cpp b/poet/src/InitGetSet.cpp deleted file mode 100644 index 1f5a813a..00000000 --- a/poet/src/InitGetSet.cpp +++ /dev/null @@ -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 -#include -#include -#include -#include -#include - -std::vector PhreeqcInit::dump_solution_names(int cell_number) { - constexpr std::size_t SKIP_H_O_CB = 3; - - std::vector 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 union_sol_names, - essential_names &names) { - // Exchange - if (this->Get_exchange(cell_number) != NULL) { - auto *exchange = this->Get_exchange(cell_number); - - std::vector 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 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 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 comp_formulas; - std::vector 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 PhreeqcInit::find_all_valence_states( - const std::vector &solution_names, const std::size_t offset) { - std::vector solution_with_valences( - solution_names.begin(), solution_names.begin() + offset); - - // to avoid duplicates store already evaluated master species - std::set 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 -PhreeqcInit::get_essential_values_init(std::size_t cell_number, - const std::vector &order) { - std::vector essentials; - - // Solutions - if (this->Get_solution(cell_number) != NULL) { - SolutionWrapper solution(this->Get_solution(cell_number), order); - std::vector sol_values(solution.size()); - std::span 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 exch_values(exchange.size()); - std::span 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 kin_values(kinetic.size()); - std::span 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 equ_values(equilibrium.size()); - std::span 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 surf_values(surface.size()); - std::span surf_span{surf_values}; - - surface.get(surf_span); - - essentials.insert(essentials.end(), surf_values.begin(), surf_values.end()); - } - - return essentials; -} - -std::map PhreeqcInit::raw_dumps() { - std::map 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; -} \ No newline at end of file diff --git a/poet/src/PhreeqcMatrix/Access.cpp b/poet/src/PhreeqcMatrix/Access.cpp new file mode 100644 index 00000000..d9f8eb61 --- /dev/null +++ b/poet/src/PhreeqcMatrix/Access.cpp @@ -0,0 +1,247 @@ +#include +#include +#include +#include +#include +#include +#include + +PhreeqcMatrix PhreeqcMatrix::subset(const std::vector &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 &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 &indices) { +// std::map> new_map; +// std::map> 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 &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(PhreeqcMatrix::PhreeqcComponent::SURFACE_COMPS); + component++) { + std::vector values; + for (const auto &[id, elements] : _m_map) { + std::vector names; + for (const auto &element : elements) { + if (static_cast(element.type) == component) { + names.push_back(element.name); + } + } + std::vector 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::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::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 +PhreeqcMatrix::getSolutionNames(bool include_h_o_charge) const { + std::vector 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 +static inline std::vector get_names_from_internal_vector( + const std::map> &map, int id) { + + const auto &it = map.find(id); + + if (it == map.end()) { + return {}; + } + + std::vector names; + + for (const auto &element : it->second) { + if (element.type == comp) { + names.push_back(element.name); + } + } + + return names; +} + +std::vector PhreeqcMatrix::getExchanger(int id) const { + return get_names_from_internal_vector( + this->_m_internal_names, id); +} + +std::vector PhreeqcMatrix::getKineticsNames(int id) const { + return get_names_from_internal_vector( + this->_m_internal_names, id); +} + +std::vector PhreeqcMatrix::getEquilibriumNames(int id) const { + return get_names_from_internal_vector( + this->_m_internal_names, id); +} + +std::vector PhreeqcMatrix::getSurfaceCompNames(int id) const { + return get_names_from_internal_vector( + this->_m_internal_names, id); +} + +std::vector PhreeqcMatrix::getSurfaceChargeNames(int id) const { + return get_names_from_internal_vector( + this->_m_internal_names, id); +} + +std::vector PhreeqcMatrix::getSolutionPrimaries() const { + return std::vector(_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; +} \ No newline at end of file diff --git a/poet/src/PhreeqcMatrix/Ctor.cpp b/poet/src/PhreeqcMatrix/Ctor.cpp new file mode 100644 index 00000000..b31b6632 --- /dev/null +++ b/poet/src/PhreeqcMatrix/Ctor.cpp @@ -0,0 +1,48 @@ +#include "IPhreeqc.hpp" +#include "PhreeqcMatrix.hpp" + +#include +#include +#include +#include + +PhreeqcMatrix::PhreeqcMatrix(const std::string &database, + const std::string &input_script) + : _m_database(database) { + this->_m_pqc = std::make_shared(); + + 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; +} \ No newline at end of file diff --git a/poet/src/PhreeqcMatrix/Init.cpp b/poet/src/PhreeqcMatrix/Init.cpp new file mode 100644 index 00000000..cf760c38 --- /dev/null +++ b/poet/src/PhreeqcMatrix/Init.cpp @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#include + +constexpr std::size_t SKIP_H_O_CB = 3; + +static std::vector dump_solution_names(cxxSolution *solution, + Phreeqc *phreeqc) { + std::vector placeholder; + + return phreeqc->find_all_valence_states( + SolutionWrapper::names(solution, placeholder), SKIP_H_O_CB); +} + +template +static void +base_add_to_element_vector(const T &wrapper, + const std::vector &names, + std::vector &elements) { + std::vector values(wrapper.size()); + std::span 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 +static void ex_ki_eq_add_to_element_vector( + T *phreeqc_component, std::vector &elements, + std::vector &b_names) { + if (phreeqc_component == NULL) { + return; + } + + std::vector ctor_names; + std::vector names = + Wrapper_t::names(phreeqc_component, ctor_names); + + for (const auto &name : ctor_names) { + b_names.push_back( + {static_cast(comp), name}); + } + + Wrapper_t wrapper(phreeqc_component, ctor_names); + + base_add_to_element_vector(wrapper, names, elements); +} + +static void +surface_add_to_element_vector(Phreeqc *phreeqc, cxxSurface *surface, + const std::vector &solution_names, + std::vector &elements, + std::vector &b_names, + std::set &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 comp_formulas; + std::vector charge_names; + + std::vector 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::SURACE_COMP), + comp}); + } + + for (const auto &charge : charge_names) { + b_names.push_back( + {static_cast( + PhreeqcMatrix::base_names::Components::SURFACE_CHARGE), + charge}); + } + + SurfaceWrapper wrapper(surface, surface_primaries, comp_formulas, + charge_names); + + base_add_to_element_vector( + wrapper, surf_names, elements); +} + +static std::pair, + std::vector> +create_vector_from_phreeqc(Phreeqc *phreeqc, int id, + const std::vector &solution_names, + std::set &surface_primaries) { + std::vector elements; + std::vector 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( + 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( + 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( + Utilities::Rxn_find(phreeqc->Get_Rxn_kinetics_map(), id), elements, + b_names); + // Equilibrium + ex_ki_eq_add_to_element_vector( + 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 find_all_solutions(Phreeqc *phreeqc) { + std::vector> all_names; + + for (auto &[id, solution] : phreeqc->Get_Rxn_solution_map()) { + if (id < 0) { + continue; + } + std::vector names = dump_solution_names(&solution, phreeqc); + all_names.push_back(names); + } + + std::vector union_names; + + for (const auto &vec : all_names) { + std::vector 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 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> elements_to_remove; + + for (const auto &[id, elements] : _m_map) { + elements_to_remove.push_back({}); + std::vector &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 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()); + } + } +} \ No newline at end of file diff --git a/poet/src/PhreeqcMatrix/Misc.cpp b/poet/src/PhreeqcMatrix/Misc.cpp new file mode 100644 index 00000000..fffe8ad8 --- /dev/null +++ b/poet/src/PhreeqcMatrix/Misc.cpp @@ -0,0 +1,53 @@ +#include "PhreeqcMatrix.hpp" +#include + +std::vector PhreeqcMatrix::getIds() const { + std::vector ids; + for (const auto &[id, _] : _m_map) { + ids.push_back(id); + } + + return ids; +} + +// std::array PhreeqcMatrix::getComponentCount(int cell_id) +// const { +// std::array 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(element.type)]++; +// } + +// return counts; +// } + +std::map PhreeqcMatrix::getDumpStringsPQI() const { + std::map 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; } \ No newline at end of file diff --git a/poet/src/Wrapper/EquilibriumCompWrapper.cpp b/poet/src/Wrapper/EquilibriumCompWrapper.cpp index 2a7c5672..e2dee0a0 100644 --- a/poet/src/Wrapper/EquilibriumCompWrapper.cpp +++ b/poet/src/Wrapper/EquilibriumCompWrapper.cpp @@ -23,7 +23,7 @@ std::vector EquilibriumWrapper::EquilibriumCompWrapper::names( std::vector 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; diff --git a/poet/test/CMakeLists.txt b/poet/test/CMakeLists.txt new file mode 100644 index 00000000..e5057fe3 --- /dev/null +++ b/poet/test/CMakeLists.txt @@ -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) \ No newline at end of file diff --git a/poet/test/barite_db.dat b/poet/test/barite_db.dat new file mode 100644 index 00000000..0d6b6fd0 --- /dev/null +++ b/poet/test/barite_db.dat @@ -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 \ No newline at end of file diff --git a/poet/test/barite_het.pqi b/poet/test/barite_het.pqi new file mode 100644 index 00000000..b8cf1e10 --- /dev/null +++ b/poet/test/barite_het.pqi @@ -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 \ No newline at end of file diff --git a/poet/test/testPhreeqcMatrix.cpp b/poet/test/testPhreeqcMatrix.cpp new file mode 100644 index 00000000..cf6b00f5 --- /dev/null +++ b/poet/test/testPhreeqcMatrix.cpp @@ -0,0 +1,204 @@ +#include +#include +#include + +#include + +#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({"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({"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({"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({"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 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 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")); +} \ No newline at end of file diff --git a/poet/test/testPhreeqcMatrix.hpp.in b/poet/test/testPhreeqcMatrix.hpp.in new file mode 100644 index 00000000..f22341d1 --- /dev/null +++ b/poet/test/testPhreeqcMatrix.hpp.in @@ -0,0 +1,105 @@ +#pragma once + +#include +#include +#include + +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 expected_names = { + "ID", "H", "O", "Charge", "Ca", "Cl", + "Mg", "Na", "Calcite_eq", "Calcite_si", "Dolomite_eq", "Dolomite_si"}; + +const std::vector 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 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 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 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::quiet_NaN(), + std::numeric_limits::quiet_NaN(), + std::numeric_limits::quiet_NaN(), + std::numeric_limits::quiet_NaN(), + 0.99937952729135193, + 0}; + +const std::vector expected_errors = { + 0, + 1e-3, + 1e-3, + 1e-15, + 1e-5, + 1e-5, + 1e-5, + 1e-5, + 1e-5, + std::numeric_limits::quiet_NaN(), + std::numeric_limits::quiet_NaN(), + std::numeric_limits::quiet_NaN(), + std::numeric_limits::quiet_NaN(), + 1e-5, + 0}; + +const std::vector expected_names_erased = { + "ID", "H", "O", "Charge", "Ba", "Cl", "S(-2)", + "S(6)", "Sr", "Barite", "Barite_p1", "Celestite", "Celestite_p1"}; + +const std::vector expected_names_subset = { + "ID", "H", "O", "Charge", "Ba", "Cl", + "S(-2)", "S(6)", "Sr", "Celestite_eq", "Celestite_si"}; +} // namespace barite_test \ No newline at end of file diff --git a/src/IPhreeqc.hpp b/src/IPhreeqc.hpp index a6daa5f1..36c3a6ae 100644 --- a/src/IPhreeqc.hpp +++ b/src/IPhreeqc.hpp @@ -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); diff --git a/src/phreeqcpp/GFZ.cpp b/src/phreeqcpp/GFZ.cpp new file mode 100644 index 00000000..6a8aee5f --- /dev/null +++ b/src/phreeqcpp/GFZ.cpp @@ -0,0 +1,75 @@ +#include "Phreeqc.h" + +std::vector Phreeqc::find_all_valence_states( + const std::vector &&solution_names, const std::size_t offset) { + std::vector solution_with_valences( + solution_names.begin(), solution_names.begin() + offset); + + // to avoid duplicates store already evaluated master species + std::set 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; +} \ No newline at end of file diff --git a/src/phreeqcpp/Phreeqc.h b/src/phreeqcpp/Phreeqc.h index 20371f8d..75b4b031 100644 --- a/src/phreeqcpp/Phreeqc.h +++ b/src/phreeqcpp/Phreeqc.h @@ -1856,14 +1856,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 keycount; // used to mark keywords that have been read + std::vector 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 + find_all_valence_states(const std::vector &&solution_names, + const std::size_t offset); + + static const class const_iso iso_defaults[]; + static const int count_iso_defaults; }; #endif /* _INC_PHREEQC_H */