Compare commits

...

5 Commits

9 changed files with 420 additions and 120 deletions

View File

@ -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

View File

@ -7,6 +7,7 @@
#include <vector>
#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<std::string> Whole vector of names. Empty if no KINETICS
* is defined
* @return std::vector<std::string> Whole vector of names. Empty if no
* KINETICS is defined
*/
std::vector<std::string> getMatrixKinetics() const;
std::vector<std::string> getMatrixKinetics() const;
/**
* @brief Returns all column names of the Matrix pertaining to EQUILIBRIUM
*
* This function returns a string vector.
*
* @return std::vector<std::string> Whole vector of names. Empty if no EQUILIBRIUM
* is defined
* @return std::vector<std::string> Whole vector of names. Empty if no
* EQUILIBRIUM is defined
*/
std::vector<std::string> getMatrixEquilibrium() const;
std::vector<std::string> getMatrixEquilibrium() const;
/*
/*
@brief Returns all column names of the Matrix pertaining to
quantities that must be transported
@brief Returns all column names of the Matrix pertaining to
quantities that must be transported
@return std::vector<std::string> vector of names
@return std::vector<std::string> vector of names
*/
std::vector<std::string> 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<std::string> vector of names
*/
std::vector<std::string> 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<std::string> vector of names
*/
std::vector<std::string> 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<std::string> vector of names
*/
std::vector<std::string> getMatrixOutOnly() const;
std::vector<std::string> getSelectedOutputNames() const {
return this->_m_selected_output_parser->getHeader();
}
private:
std::map<int, std::vector<element>> _m_map;
@ -367,6 +384,7 @@ private:
std::shared_ptr<IPhreeqc> _m_pqc;
std::shared_ptr<PhreeqcKnobs> _m_knobs;
std::shared_ptr<PhreeqcSelectedOutputParser> _m_selected_output_parser;
std::string _m_database;

View File

@ -0,0 +1,95 @@
#include "IPhreeqc.hpp"
#include <cstdint>
/**
* @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<std::string> 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<double> A vector of output values corresponding to the
* specified cell.
*/
std::vector<double> 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<std::string> _m_headings;
};

View File

@ -1,6 +1,7 @@
#include "PhreeqcEngine.hpp"
#include <cstddef>
#include <iomanip>
#include <memory>
#include <regex>
#include <span>
#include <sstream>
@ -52,11 +53,14 @@ public:
std::unique_ptr<KineticWrapper> kineticsWrapperPtr;
std::unique_ptr<EquilibriumWrapper> equilibriumWrapperPtr;
std::unique_ptr<SurfaceWrapper> surfaceWrapperPtr;
std::unique_ptr<PhreeqcSelectedOutputParser> _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<std::string> 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<PhreeqcSelectedOutputParser>(
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<double> &data) {
@ -209,6 +220,20 @@ void PhreeqcEngine::Impl::get_essential_values(std::span<double> &data) {
std::span<double> surf_span{
data.subspan(offset, this->surfaceWrapperPtr->size())};
this->surfaceWrapperPtr->get(surf_span);
offset += this->surfaceWrapperPtr->size();
}
if (this->has_selected_output) {
std::vector<double> selected_output_values =
this->_m_selected_output_parser->getValues(1);
std::span<double> 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();
}
}

View File

@ -85,8 +85,8 @@ PhreeqcMatrix::STLExport PhreeqcMatrix::get(VectorExportType type,
}
for (std::size_t component = 1;
component <=
static_cast<std::size_t>(PhreeqcMatrix::PhreeqcComponent::SURFACE_COMPS);
component <= static_cast<std::size_t>(
PhreeqcMatrix::PhreeqcComponent::SELECTED_OUTPUT);
component++) {
std::vector<std::string> 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<std::string> PhreeqcMatrix::getMatrixKinetics() const {
std::vector<std::string> names;
std::vector<std::string> 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);
}
}
}
}
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;
}
return names;
}
// MDL
std::vector<std::string> PhreeqcMatrix::getMatrixEquilibrium() const {
std::vector<std::string> names;
std::vector<std::string> 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);
}
}
}
}
std::vector<std::string> names;
std::vector<std::string> 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<std::string> PhreeqcMatrix::getMatrixTransported() const {
std::vector<std::string> names;
std::vector<std::string> names;
const std::vector<std::string> to_remove = {
"tc", "patm", "SolVol", "pH", "pe"
};
const std::vector<std::string> 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);
}
// 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<std::string> PhreeqcMatrix::getMatrixOutOnly() const {
// MDL we must append here selected_output / user_punch
std::vector<std::string> defaultnames = {
"tc", "patm", "SolVol", "pH", "pe"
};
std::vector<std::string> ret;
for (auto nm : defaultnames) {
ret.push_back(nm);
}
return ret;
// MDL we must append here selected_output / user_punch
std::vector<std::string> defaultnames = {"tc", "patm", "SolVol", "pH", "pe"};
std::vector<std::string> ret;
for (auto nm : defaultnames) {
ret.push_back(nm);
}
return ret;
}

View File

@ -4,7 +4,6 @@
#include <Phreeqc.h>
#include <Solution.h>
#include <cmath>
#include <memory>
#include <string>
@ -16,6 +15,9 @@ PhreeqcMatrix::PhreeqcMatrix(const std::string &database,
this->_m_pqc = std::make_shared<IPhreeqc>();
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<PhreeqcSelectedOutputParser>(this->_m_pqc.get(),
input_script);
this->_m_knobs =
std::make_shared<PhreeqcKnobs>(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;
// }

View File

@ -12,6 +12,7 @@
#include <cmath>
#include <cstddef>
#include <iterator>
#include <limits>
#include <map>
#include <set>
#include <string>
@ -224,9 +225,23 @@ 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<double> selected_output_values =
this->_m_selected_output_parser->getValues(id);
const auto &headers = this->_m_selected_output_parser->getHeader();
for (const auto &hdr_str : headers) {
elements.push_back(
{hdr_str, PhreeqcMatrix::PhreeqcComponent::SELECTED_OUTPUT, 0});
base_names.push_back(
{PhreeqcMatrix::base_names::Components::SELECTED_OUTPUT, hdr_str});
}
}
_m_map[id] = elements;
_m_internal_names[id] = base_names;
}

View File

@ -0,0 +1,170 @@
#include "PhreeqcSelectedOutputParser.hpp"
#include <limits>
#include <regex>
#include <sstream>
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();
if (this->_m_headings.size() != this->getValues(1).size()) {
throw std::runtime_error(
"Number of headings does not match number of values in selected "
"output.");
}
}
}
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<double>
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<double> 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<double>::quiet_NaN()); // If conversion fails,
// store NaN
}
}
return values;
}

View File

@ -1,4 +1,4 @@
// Time-stamp: "Last modified 2025-07-28 20:14:08 delucia"
// Time-stamp: "Last modified 2025-08-01 10:54:30 delucia"
#include <iostream>
#include <iomanip>
#include <linux/limits.h>
@ -96,6 +96,10 @@ int main(int argc, char *argv[]) {
auto outonly = pqc_mat1.getMatrixOutOnly();
std::cout << ":: pqc_mat1.getMatrixOutOnly(): \n";
std::cout << outonly << "\n\n";
auto selout = pqc_mat1.getSelectedOutputNames();
std::cout << ":: pqc_mat1.getSelectedOutputNames(): \n";
std::cout << selout << "\n\n";
return 0;
}