From 1b12899ffb4ab927c050ff3fdae865af2ed34e56 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Wed, 30 Jul 2025 15:45:33 +0200 Subject: [PATCH] Selected Output: First prototype with successfull exisiting tests --- poet/CMakeLists.txt | 1 + poet/include/PhreeqcMatrix.hpp | 72 +++++--- poet/include/PhreeqcSelectedOutputParser.hpp | 95 +++++++++++ poet/src/Engine.cpp | 25 +++ poet/src/PhreeqcMatrix/Access.cpp | 114 +++++++------ poet/src/PhreeqcMatrix/Ctor.cpp | 38 +---- poet/src/PhreeqcMatrix/Init.cpp | 20 ++- poet/src/SelectedOutputParser.cpp | 165 +++++++++++++++++++ 8 files changed, 411 insertions(+), 119 deletions(-) create mode 100644 poet/include/PhreeqcSelectedOutputParser.hpp create mode 100644 poet/src/SelectedOutputParser.cpp diff --git a/poet/CMakeLists.txt b/poet/CMakeLists.txt index 5050dca7..0103266a 100644 --- a/poet/CMakeLists.txt +++ b/poet/CMakeLists.txt @@ -2,6 +2,7 @@ set(POET_SOURCE_FILES src/Engine.cpp src/Runner.cpp src/Knobs.cpp + src/SelectedOutputParser.cpp #Wrappers src/Wrapper/EquilibriumWrapper.cpp src/Wrapper/EquilibriumCompWrapper.cpp diff --git a/poet/include/PhreeqcMatrix.hpp b/poet/include/PhreeqcMatrix.hpp index d95b71a7..765cbd1b 100644 --- a/poet/include/PhreeqcMatrix.hpp +++ b/poet/include/PhreeqcMatrix.hpp @@ -7,6 +7,7 @@ #include #include "PhreeqcKnobs.hpp" +#include "PhreeqcSelectedOutputParser.hpp" class IPhreeqc; @@ -170,7 +171,8 @@ public: EXCHANGE, KINETIC, EQUILIBRIUM, - SURFACE_COMPS + SURFACE_COMPS, + SELECTED_OUTPUT }; struct element { @@ -185,12 +187,17 @@ public: KINETICS, EQUILIBRIUM, SURACE_COMP, - SURFACE_CHARGE + SURFACE_CHARGE, + SELECTED_OUTPUT } type; std::string name; }; + const PhreeqcSelectedOutputParser &getSelectedOutput() const { + return *_m_selected_output_parser; + } + /** * @brief Get all found solution names. * @@ -319,41 +326,51 @@ public: * * This function returns a string vector. * - * @return std::vector Whole vector of names. Empty if no KINETICS - * is defined + * @return std::vector Whole vector of names. Empty if no + * KINETICS is defined */ - std::vector getMatrixKinetics() const; - + std::vector getMatrixKinetics() const; + /** * @brief Returns all column names of the Matrix pertaining to EQUILIBRIUM * * This function returns a string vector. * - * @return std::vector Whole vector of names. Empty if no EQUILIBRIUM - * is defined + * @return std::vector Whole vector of names. Empty if no + * EQUILIBRIUM is defined */ - std::vector getMatrixEquilibrium() const; + std::vector getMatrixEquilibrium() const; - - /* - - @brief Returns all column names of the Matrix pertaining to - quantities that must be transported - - @return std::vector vector of names + /* + @brief Returns all column names of the Matrix pertaining to + quantities that must be transported + + @return std::vector vector of names + + */ + std::vector getMatrixTransported() const; + + /* + + @brief Returns all column names of the Matrix pertaining to + quantities that must NOT be transported but have to be included in + the output + + @return std::vector vector of names + */ + std::vector getMatrixOutOnly() const; + + /** + * @brief Returns all column names of the Matrix pertaining to + * quantities that are user defined in the SELECTED_OUTPUT or USER_PUNCH + * blocks of the Phreeqc input script. + * + * @return std::vector vector of names */ - std::vector getMatrixTransported() const; - - /* - - @brief Returns all column names of the Matrix pertaining to - quantities that must NOT be transported but have to be included in - the output - - @return std::vector vector of names - */ - std::vector getMatrixOutOnly() const; + std::vector getSelectedOutputNames() const { + return this->_m_selected_output_parser->getHeader(); + } private: std::map> _m_map; @@ -367,6 +384,7 @@ private: std::shared_ptr _m_pqc; std::shared_ptr _m_knobs; + std::shared_ptr _m_selected_output_parser; std::string _m_database; diff --git a/poet/include/PhreeqcSelectedOutputParser.hpp b/poet/include/PhreeqcSelectedOutputParser.hpp new file mode 100644 index 00000000..1d00d4c3 --- /dev/null +++ b/poet/include/PhreeqcSelectedOutputParser.hpp @@ -0,0 +1,95 @@ +#include "IPhreeqc.hpp" +#include + +/** + * @class PhreeqcSelectedOutputParser + * @brief Parses and manages the SELECTED_OUTPUT block from a PHREEQC input + * script. + * + * This class is responsible for extracting and handling the SELECTED_OUTPUT + * block from a PHREEQC input script. It provides methods to check if a + * SELECTED_OUTPUT block is present, retrieve its string representation, access + * the parsed header, and obtain output values for a specific cell. + * + * @note The class requires an instance of IPhreeqc to function. + */ +class PhreeqcSelectedOutputParser { +public: + /** + * @brief Constructs a PhreeqcSelectedOutputParser object. + * + * Initializes the parser with a given IPhreeqc instance and an input script. + * The parser will use the provided PHREEQC instance to process the specified + * input script and extract selected output data. + * + * @param pqc_instance Pointer to an IPhreeqc instance used for PHREEQC + * calculations. + * @param input_script The PHREEQC input script to be processed. + */ + PhreeqcSelectedOutputParser(IPhreeqc *pqc_instance, + const std::string &input_script); + + /** + * @brief Constructs a PhreeqcSelectedOutputParser by associating it with an + * existing IPhreeqc instance and initializing it using another + * PhreeqcSelectedOutputParser. + * + * This constructor allows creating a new parser that shares the IPhreeqc + * context and copies relevant state or configuration from an existing parser + * instance. + * + * @param pqc_instance Pointer to an IPhreeqc instance to be used by the + * parser. + * @param input_parser Reference to an existing PhreeqcSelectedOutputParser to + * initialize from. + */ + PhreeqcSelectedOutputParser(IPhreeqc *pqc_instance, + const PhreeqcSelectedOutputParser &input_parser); + + /** + * @brief Checks if selected output is available. + * + * @return true if selected output is present, false otherwise. + */ + bool hasSelectedOutput() const { return _m_has_selected_output; } + + /** + * @brief Retrieves the string representation of the selected output block. + * + * @return A constant reference to the string containing the selected output + * block. + */ + const std::string &getSelectedOutputBlockString() const { + return _m_selected_output_block_string; + } + + /** + * @brief Retrieves the header row of the selected output. + * + * @return A vector of strings containing the column headings. + */ + std::vector getHeader() const { return _m_headings; } + + /** + * @brief Retrieves the selected output values for a specified cell. + * + * This function returns a vector containing the output values associated with + * the given cell ID. + * + * @param cell_id The identifier of the Phreeqc cell for which to retrieve the + * output values. + * @return std::vector A vector of output values corresponding to the + * specified cell. + */ + std::vector getValues(std::uint32_t cell_id) const; + +private: + void parseSelectedOutputBlock(const std::string &input_script); + void parseHeader(); + + bool _m_has_selected_output{false}; // true if selected output was defined + std::string _m_selected_output_block_string; + + IPhreeqc *_m_pqc_instance; + std::vector _m_headings; +}; diff --git a/poet/src/Engine.cpp b/poet/src/Engine.cpp index 7c81b3c3..88a41897 100644 --- a/poet/src/Engine.cpp +++ b/poet/src/Engine.cpp @@ -1,6 +1,7 @@ #include "PhreeqcEngine.hpp" #include #include +#include #include #include #include @@ -52,11 +53,14 @@ public: std::unique_ptr kineticsWrapperPtr; std::unique_ptr equilibriumWrapperPtr; std::unique_ptr surfaceWrapperPtr; + std::unique_ptr _m_selected_output_parser; bool has_exchange = false; bool has_kinetics = false; bool has_equilibrium = false; bool has_surface = false; + bool has_selected_output = false; + struct InitCell { std::vector solutions; bool with_redox; @@ -94,6 +98,10 @@ PhreeqcEngine::Impl::Impl(const PhreeqcMatrix &pqc_mat, const int cell_id) { const std::string pqc_string = replaceRawKeywordID(pqc_mat.getDumpStringsPQI(cell_id)); + this->_m_selected_output_parser = + std::make_unique( + this, pqc_mat.getSelectedOutput()); + this->RunString(pqc_string.c_str()); InitCell cell = {pqc_mat.getSolutionNames(), @@ -173,6 +181,9 @@ void PhreeqcEngine::Impl::init_wrappers(const InitCell &cell) { this->has_surface = true; } + + this->has_selected_output = + this->_m_selected_output_parser->hasSelectedOutput(); } void PhreeqcEngine::Impl::get_essential_values(std::span &data) { @@ -209,6 +220,20 @@ void PhreeqcEngine::Impl::get_essential_values(std::span &data) { std::span surf_span{ data.subspan(offset, this->surfaceWrapperPtr->size())}; this->surfaceWrapperPtr->get(surf_span); + + offset += this->surfaceWrapperPtr->size(); + } + + if (this->has_selected_output) { + std::vector selected_output_values = + this->_m_selected_output_parser->getValues(1); + + std::span sel_out_span{ + data.subspan(offset, selected_output_values.size())}; + std::copy(selected_output_values.begin(), selected_output_values.end(), + sel_out_span.begin()); + + offset += selected_output_values.size(); } } diff --git a/poet/src/PhreeqcMatrix/Access.cpp b/poet/src/PhreeqcMatrix/Access.cpp index 01d77b49..b5364119 100644 --- a/poet/src/PhreeqcMatrix/Access.cpp +++ b/poet/src/PhreeqcMatrix/Access.cpp @@ -85,8 +85,8 @@ PhreeqcMatrix::STLExport PhreeqcMatrix::get(VectorExportType type, } for (std::size_t component = 1; - component <= - static_cast(PhreeqcMatrix::PhreeqcComponent::SURFACE_COMPS); + component <= static_cast( + PhreeqcMatrix::PhreeqcComponent::SELECTED_OUTPUT); component++) { std::vector values; for (const auto &[id, elements] : _m_map) { @@ -241,79 +241,75 @@ double PhreeqcMatrix::operator()(int cell_id, const std::string &name) const { return it->value; } - // MDL std::vector PhreeqcMatrix::getMatrixKinetics() const { - std::vector names; - - auto n = this->getIds().size(); - for (auto i = 0; igetKineticsNames(i); - for (auto nam : pqc_kinnames ) { - for (auto mat_name : this->get().names){ - if (mat_name.starts_with(nam)) { - // check if we already have this mat_name - if (std::find(names.begin(), names.end(), mat_name) == names.end()) { - names.push_back(mat_name); - } - } - } - } - } - return names; -} + std::vector names; + auto n = this->getIds().size(); + for (auto i = 0; i < n; ++i) { + auto pqc_kinnames = this->getKineticsNames(i); + for (auto nam : pqc_kinnames) { + for (auto mat_name : this->get().names) { + if (mat_name.starts_with(nam)) { + // check if we already have this mat_name + if (std::find(names.begin(), names.end(), mat_name) == names.end()) { + names.push_back(mat_name); + } + } + } + } + } + return names; +} // MDL std::vector PhreeqcMatrix::getMatrixEquilibrium() const { - std::vector names; - std::vector mat_names = this->get().names; - auto n = this->getIds().size(); - for (auto i = 0; igetEquilibriumNames(i); - for (auto nam : pqc_eqnames ) { - for (auto mat_name : mat_names){ - if (mat_name.starts_with(nam)) { - // check if we already have this mat_name - if (std::find(names.begin(), names.end(), mat_name) == names.end()) { - names.push_back(mat_name); - } - } - } - } + std::vector names; + std::vector mat_names = this->get().names; + auto n = this->getIds().size(); + for (auto i = 0; i < n; ++i) { + auto pqc_eqnames = this->getEquilibriumNames(i); + for (auto nam : pqc_eqnames) { + for (auto mat_name : mat_names) { + if (mat_name.starts_with(nam)) { + // check if we already have this mat_name + if (std::find(names.begin(), names.end(), mat_name) == names.end()) { + names.push_back(mat_name); + } + } + } } - return names; + } + return names; } // MDL std::vector PhreeqcMatrix::getMatrixTransported() const { - std::vector names; - - const std::vector to_remove = { - "tc", "patm", "SolVol", "pH", "pe" - }; + std::vector names; - // sols contains all solutes; we must remove { tc, patm, SolVol, pH, pe } - auto sols = this->getSolutionNames(); - for (auto name : sols) { - if (std::find(to_remove.begin(), to_remove.end(), name) == to_remove.end()) { - names.push_back(name); - } + const std::vector to_remove = {"tc", "patm", "SolVol", "pH", + "pe"}; + + // sols contains all solutes; we must remove { tc, patm, SolVol, pH, pe } + auto sols = this->getSolutionNames(); + for (auto name : sols) { + if (std::find(to_remove.begin(), to_remove.end(), name) == + to_remove.end()) { + names.push_back(name); } - - return names; + } + + return names; } // MDL std::vector PhreeqcMatrix::getMatrixOutOnly() const { - // MDL we must append here selected_output / user_punch - std::vector defaultnames = { - "tc", "patm", "SolVol", "pH", "pe" - }; - std::vector ret; - for (auto nm : defaultnames) { - ret.push_back(nm); - } - return ret; + // MDL we must append here selected_output / user_punch + std::vector defaultnames = {"tc", "patm", "SolVol", "pH", "pe"}; + std::vector ret; + for (auto nm : defaultnames) { + ret.push_back(nm); + } + return ret; } diff --git a/poet/src/PhreeqcMatrix/Ctor.cpp b/poet/src/PhreeqcMatrix/Ctor.cpp index 0c0b4ad0..c8061167 100644 --- a/poet/src/PhreeqcMatrix/Ctor.cpp +++ b/poet/src/PhreeqcMatrix/Ctor.cpp @@ -4,7 +4,6 @@ #include #include -#include #include #include @@ -16,6 +15,9 @@ PhreeqcMatrix::PhreeqcMatrix(const std::string &database, this->_m_pqc = std::make_shared(); this->_m_pqc->LoadDatabaseString(database.c_str()); + + this->_m_pqc->SetSelectedOutputStringOn(true); + this->_m_pqc->RunString(input_script.c_str()); if (this->_m_pqc->GetErrorStringLineCount() > 0) { @@ -24,38 +26,12 @@ PhreeqcMatrix::PhreeqcMatrix(const std::string &database, throw std::runtime_error("Phreeqc script error"); } + this->_m_selected_output_parser = + std::make_shared(this->_m_pqc.get(), + input_script); + this->_m_knobs = std::make_shared(this->_m_pqc.get()->GetPhreeqcPtr()); 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; -// } diff --git a/poet/src/PhreeqcMatrix/Init.cpp b/poet/src/PhreeqcMatrix/Init.cpp index 5fdc9bd7..074fd2b4 100644 --- a/poet/src/PhreeqcMatrix/Init.cpp +++ b/poet/src/PhreeqcMatrix/Init.cpp @@ -224,9 +224,25 @@ void PhreeqcMatrix::initialize() { if (id < 0) { continue; } - const auto &[elements, base_names] = create_vector_from_phreeqc( + auto [elements, base_names] = create_vector_from_phreeqc( phreeqc, id, solutions, this->_m_surface_primaries); + if (this->_m_selected_output_parser->hasSelectedOutput()) { + // Add selected output elements + std::vector selected_output_values = + this->_m_selected_output_parser->getValues(id); + + const auto &headers = this->_m_selected_output_parser->getHeader(); + + for (std::size_t i = 0; i < selected_output_values.size(); i++) { + elements.push_back({headers[i], + PhreeqcMatrix::PhreeqcComponent::SELECTED_OUTPUT, + selected_output_values[i]}); + base_names.push_back( + {PhreeqcMatrix::base_names::Components::SELECTED_OUTPUT, + headers[i]}); + } + } _m_map[id] = elements; _m_internal_names[id] = base_names; } @@ -262,4 +278,4 @@ void PhreeqcMatrix::remove_NaNs() { elements.end()); } } -} \ No newline at end of file +} diff --git a/poet/src/SelectedOutputParser.cpp b/poet/src/SelectedOutputParser.cpp new file mode 100644 index 00000000..4e288d8b --- /dev/null +++ b/poet/src/SelectedOutputParser.cpp @@ -0,0 +1,165 @@ +#include "PhreeqcSelectedOutputParser.hpp" + +#include +#include +#include + +PhreeqcSelectedOutputParser::PhreeqcSelectedOutputParser( + IPhreeqc *pqc_instance, const std::string &input_script) + : _m_pqc_instance(pqc_instance) { + if (!_m_pqc_instance->GetSelectedOutputStringOn()) { + throw std::runtime_error("Selected output string is not enabled."); + } + + this->parseSelectedOutputBlock(input_script); + + if (this->_m_has_selected_output) { + parseHeader(); + } +} + +PhreeqcSelectedOutputParser::PhreeqcSelectedOutputParser( + IPhreeqc *pqc_instance, const PhreeqcSelectedOutputParser &input_parser) + : _m_pqc_instance(pqc_instance), + _m_has_selected_output(input_parser._m_has_selected_output), + _m_selected_output_block_string( + input_parser._m_selected_output_block_string), + _m_headings(input_parser._m_headings) { + if (_m_has_selected_output) { + this->_m_pqc_instance->SetSelectedOutputStringOn(true); + this->_m_pqc_instance->RunString( + this->_m_selected_output_block_string.c_str()); + + if (this->_m_pqc_instance->GetErrorStringLineCount() > 0) { + throw std::runtime_error("Error running selected output block!"); + } + } +} + +std::string getBlockByKeyword(const std::string &input_script, + const std::string &&keyword) { + const std::regex keyword_regex("^" + keyword + "*$"); + const std::regex block_regex(R"(^[A-Z]+.*$)"); + + bool block_found = false; + std::size_t block_start = 0; + std::size_t block_end = 0; + std::size_t current_pos = 0; + + std::istringstream input_stream(input_script); + + for (std::string line; std::getline(input_stream, line); + current_pos += line.length() + 1) { + + // remove trailing whitespaces + std::size_t first_char_pos = line.find_first_not_of(" \t\r"); + + if (first_char_pos == std::string::npos) { + continue; // Skip empty lines + } + + if (std::regex_search(line, keyword_regex)) { + block_start = current_pos; + block_found = true; + continue; + } + + if (!block_found) { + continue; + } + + if (!std::regex_search( + line.substr(first_char_pos, line.length() - first_char_pos), + block_regex)) { + continue; // Skip lines that do not start with a capitalized keyword + } + + block_end = current_pos - 1; // End of KEYWORD block + break; + } + + if (!block_found) { + return {""}; + } + + return std::string(input_script, block_start, block_end - block_start + 1); +} + +void PhreeqcSelectedOutputParser::parseSelectedOutputBlock( + const std::string &input_script) { + std::string selected_output_block = + getBlockByKeyword(input_script, "SELECTED_OUTPUT"); + + std::string user_punch_block = getBlockByKeyword(input_script, "USER_PUNCH"); + + if (selected_output_block.empty() && user_punch_block.empty()) { + return; + } + + if (selected_output_block.empty() && !user_punch_block.empty()) { + throw std::runtime_error( + "USER_PUNCH block found without a SELECTED_OUTPUT block."); + } + + this->_m_has_selected_output = true; + this->_m_selected_output_block_string = + selected_output_block + user_punch_block; +} + +void PhreeqcSelectedOutputParser::parseHeader() { + const std::string selected_output_string = + _m_pqc_instance->GetSelectedOutputString(); + + std::istringstream stream(selected_output_string); + std::string header_line; + + if (!std::getline(stream, header_line)) { + throw std::runtime_error("No headings found in selected output string."); + } + + std::istringstream header_stream(header_line); + + for (std::string heading; std::getline(header_stream, heading, '\t');) { + std::size_t first_char_pos = heading.find_first_not_of(" "); + _m_headings.push_back( + heading.substr(first_char_pos, heading.length() - first_char_pos) + + "_SO"); + } +} + +std::vector +PhreeqcSelectedOutputParser::getValues(std::uint32_t cell_id) const { + if (!this->hasSelectedOutput()) { + return {}; // No selected output defined + } + // this->_m_pqc_instance->SetCurrentSelectedOutputUserNumber(cell_id); + + const std::string selected_output_string = + _m_pqc_instance->GetSelectedOutputString(); + + if (selected_output_string.empty()) { + throw std::runtime_error("Selected output string is empty."); + } + + std::istringstream stream(selected_output_string); + std::string line, tmp_line; + std::vector values; + + while (std::getline(stream, tmp_line)) { + line = tmp_line; + } + + std::istringstream line_stream(line); + + for (std::string value; std::getline(line_stream, value, '\t');) { + try { + values.push_back(std::stod(value)); + } catch (const std::invalid_argument &) { + values.push_back( + std::numeric_limits::quiet_NaN()); // If conversion fails, + // store NaN + } + } + + return values; +}