Merge branch 'golem' into 'poet'

Merge features from golem

See merge request naaice/iphreeqc!30

poet/CMakeLists.txt

@@ -61,5 +61,8 @@ if (BUILD_TESTING AND STANDALONE_BUILD)
     gtest_discover_tests(poet_test)
 endif()
 
-add_executable(golemrunner test/testGolemRunner.cpp)
-target_link_libraries(golemrunner IPhreeqcPOET)
+add_executable(testGolemRunner test/testGolemRunner.cpp)
+target_link_libraries(testGolemRunner IPhreeqcPOET)
+
+add_executable(testGetters     test/testGetters.cpp)
+target_link_libraries(testGetters IPhreeqcPOET)

poet/include/PhreeqcMatrix.hpp

@@ -313,6 +313,48 @@ public:
    */
   bool withRedox() const { return _m_with_redox; }
 
+  // MDL
+  /**
+   * @brief Returns all column names of the Matrix pertaining to KINETICS
+   *
+   * This function returns a string vector.
+   *
+   * @return std::vector<std::string> Whole vector of names. Empty if no KINETICS
+   * is defined
+   */
+   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
+   */
+   std::vector<std::string> getMatrixEquilibrium() const;
+
+      
+   /*
+   
+   @brief Returns all column names of the Matrix pertaining to
+   quantities that must be transported
+   
+   @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;
+
 private:
   std::map<int, std::vector<element>> _m_map;
   std::map<int, std::vector<base_names>> _m_internal_names;
@@ -330,4 +372,4 @@ private:
 
   bool _m_with_h0_o0;
   bool _m_with_redox;
-};
+};

poet/src/Engine.cpp

@@ -130,6 +130,12 @@ void PhreeqcEngine::Impl::run(double time_step) {
   const std::string runs_string =
       "RUN_CELLS\n -cells 1\n -time_step " + time_ss.str() + "\nEND\n";
   this->RunString(runs_string.c_str());
+
+  if (this->GetErrorStringLineCount() > 0) {
+    std::cerr << ":: Error in Phreeqc script: " << this->GetErrorString()
+              << "\n";
+    throw std::runtime_error("Phreeqc script error");
+  }
 }
 
 void PhreeqcEngine::Impl::init_wrappers(const InitCell &cell) {
@@ -242,4 +248,4 @@ void PhreeqcEngine::Impl::set_essential_values(const std::span<double> &data) {
         data.subspan(offset, this->surfaceWrapperPtr->size())};
     this->surfaceWrapperPtr->set(surf_span);
   }
-}
+}

poet/src/PhreeqcMatrix/Access.cpp

@@ -119,6 +119,7 @@ PhreeqcMatrix::STLExport PhreeqcMatrix::get(VectorExportType type,
     for (; column_index < result.names.size(); column_index++) {
       for (const auto &[_, elements] : _m_map) {
         double value_to_add = std::numeric_limits<double>::quiet_NaN();
+        // double value_to_add;
         for (const auto &curr_element : elements) {
           const std::string &curr_element_name = curr_element.name;
 
@@ -238,4 +239,81 @@ 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;
+    
+    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<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);
+		    }
+		}
+	    }
+	}
+    }
+    return names;
+}
+
+// MDL
+std::vector<std::string> PhreeqcMatrix::getMatrixTransported() const {
+    std::vector<std::string> names;
+    
+    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);
+	}
+    }
+    
+    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;
+}

poet/src/PhreeqcMatrix/Ctor.cpp

@@ -18,6 +18,12 @@ PhreeqcMatrix::PhreeqcMatrix(const std::string &database,
   this->_m_pqc->LoadDatabaseString(database.c_str());
   this->_m_pqc->RunString(input_script.c_str());
 
+  if (this->_m_pqc->GetErrorStringLineCount() > 0) {
+    std::cerr << ":: Error in Phreeqc script: "
+              << this->_m_pqc->GetErrorString() << "\n";
+    throw std::runtime_error("Phreeqc script error");
+  }
+
   this->_m_knobs =
       std::make_shared<PhreeqcKnobs>(this->_m_pqc.get()->GetPhreeqcPtr());
 
@@ -52,4 +58,4 @@ PhreeqcMatrix::PhreeqcMatrix(const std::string &database,
 //   _m_database = other._m_database;
 
 //   return *this;
-// }
+// }

poet/src/Wrapper/KineticCompWrapper.cpp

@@ -30,11 +30,11 @@ KineticWrapper::KineticCompWrapper::names(const cxxKineticsComp &comp) {
   std::vector<std::string> names;
 
   const std::string &comp_name = comp.Get_rate_name();
-  names.push_back(comp_name);
+  names.push_back(comp_name + "_kin");
 
   for (std::size_t i = 0; i < comp.Get_d_params().size(); i++) {
     names.push_back(comp_name + "_p" + std::to_string(i + 1));
   }
 
   return names;
-}
+}

poet/src/Wrapper/SolutionWrapper.cpp

@@ -18,6 +18,11 @@ void SolutionWrapper::get(std::span<LDBLE> &data) const {
   data[0] = solution->Get_total_h();
   data[1] = solution->Get_total_o();
   data[2] = solution->Get_cb();
+  data[3] = solution->Get_tc();
+  data[4] = solution->Get_patm();
+  data[5] = solution->Get_soln_vol();
+  data[6] = solution->Get_ph();
+  data[7] = solution->Get_pe();
 
   const cxxNameDouble &totals =
       (_with_redox ? solution->Get_totals()
@@ -41,6 +46,8 @@ void SolutionWrapper::set(const std::span<LDBLE> &data) {
   const double &total_h = data[0];
   const double &total_o = data[1];
   const double &cb = data[2];
+  const double &tc = data[3];
+  const double &patm = data[4];
 
   for (const auto &tot_name : solution_order) {
     const double value = data[i++];
@@ -51,7 +58,7 @@ void SolutionWrapper::set(const std::span<LDBLE> &data) {
     new_totals[tot_name] = value;
   }
 
-  this->solution->Update(total_h, total_o, cb,
+  this->solution->Update(total_h, total_o, cb, tc, patm,
                          _with_redox ? new_totals
                                      : new_totals.Simplify_redox());
 }
@@ -85,4 +92,4 @@ SolutionWrapper::names(cxxSolution *solution, bool include_h0_o0,
 
   names.insert(names.end(), names_set.begin(), names_set.end());
   return names;
-}
+}

poet/src/Wrapper/SolutionWrapper.hpp

@@ -27,10 +27,12 @@ private:
   cxxSolution *solution;
   const std::vector<std::string> solution_order;
 
-  static constexpr std::array<const char *, 3> ESSENTIALS = {"H", "O",
-                                                             "Charge"};
+  static constexpr std::array<const char *, 8> ESSENTIALS = {
+      "H",      "O",  "Charge", "tc", "patm",
+
+      "SolVol", "pH", "pe"}; // MDL; ML: only output
 
   static constexpr std::size_t NUM_ESSENTIALS = ESSENTIALS.size();
 
   const bool _with_redox;
-};
+};

poet/test/run_kin_cor_end2.pqi

@@ -0,0 +1,110 @@
+SOLUTION  1
+units	mol/kgw
+temp	35
+pH      7
+pe      4 
+Cl      0.003 charge
+Na      0.003
+PURE  1
+Calcite 0.0  10
+
+SOLUTION 2
+units   mol/kgw
+pH     6.921 
+pe    -1.258
+water  1
+temp  17
+Na   1e-15
+Cl   1e-15 charge
+Ca   1e-15
+C    1e-15
+Fe   1e-15
+S    1e-15
+Al   1e-15
+Si   7.963e-5
+
+PURE  2
+Pyrite  0.0  2.1e-4
+
+KINETICS  2
+Quartz
+  -m         2.2e-4
+  -parms     0.01  1
+
+SOLUTION 3
+units   mol/kgw
+pH     9.963
+pe    -6.769 
+temp 17
+Na   1e-15
+Cl   1e-15 
+Ca   1.345e-4
+C    1.345e-4
+Fe   1.433e-8
+S    2.866e-8
+Al   2.692e-5
+Si   2.692e-5
+
+PURE  3
+Pyrite         0 2e-4
+Calcite        0 1e-4
+
+KINETICS  3
+Kaolinite
+  -m          1.4e-3
+  -parms      0.6 1000
+Siderite
+  -m          1.2e-12
+  -parms      1.2 10
+
+SOLUTION  4 
+units	mol/kgw
+temp	17
+pH      5.576
+pe      4      O2(g) -1
+Cl      0.0003 charge
+Na      0.0003
+C       1e-15 CO2(g) -3.38
+Ca	1e-15
+Fe      1e-15
+S       1e-15
+Al      1e-15
+Si      1e-15
+
+
+SOLUTION 5 
+ units   mol/kgw
+ pH    6.653
+ pe    -3.167
+ water  1
+ temp  17
+ Na   1e-15
+ Cl   1e-15 charge
+ Ca   1e-15
+ C    1e-15
+ Fe   4.386e-9
+ S    8.772e-9
+ Al   5.031e-7
+ Si   5.031e-7
+
+PURE  5
+Pyrite  0.0  2.5e-4
+
+KINETICS 5
+Kaolinite
+  -m       1.5e-4 
+  -parms   0.6 100000
+
+SOLUTION  6
+units	mol/kgw
+temp	35
+pH      7
+pe      4 
+Cl      0.003 charge
+Na      0.003
+PURE  6
+Calcite 0.0  10
+
+RUN_CELLS
+ -cells 1 2 3 4 5 6
+END

poet/test/testGetters.cpp

@@ -0,0 +1,104 @@
+//  Time-stamp: "Last modified 2025-07-28 20:14:08 delucia"
+#include <iostream>
+#include <iomanip>
+#include <linux/limits.h>
+#include <sstream>
+#include <fstream>
+#include <algorithm>
+#include <iterator>
+#include <vector>
+#include <string>
+#include <memory>
+#include <cmath>
+
+#include "PhreeqcMatrix.hpp"
+#include "PhreeqcEngine.hpp"
+#include "PhreeqcRunner.hpp"
+
+
+std::string readFile(const std::string &path) {
+  std::string string_rpath(PATH_MAX, '\0');
+
+  if (realpath(path.c_str(), string_rpath.data()) == nullptr) {
+    throw std::runtime_error(":: Failed to resolve the realpath to file " + path);
+  }
+
+  std::ifstream file(string_rpath);
+
+  if (!file.is_open()) {
+    throw std::runtime_error(":: Failed to open file: " + path);
+  }
+
+  std::stringstream buffer;
+  buffer << file.rdbuf();
+
+  return buffer.str();
+}
+
+// pretty print a vector, standard implementation from stackoverflow 
+template<typename T>
+std::ostream & operator<<(std::ostream & os, std::vector<T> vec)
+{
+    os << "{ ";
+    std::copy(vec.begin(), vec.end(), std::ostream_iterator<T>(os, ", "));
+    os << " }";
+    return os;
+}
+
+
+int main(int argc, char *argv[]) {
+    
+    if (argc < 3) {
+	std::cout << "::" << argv[0] << ": two args needed, script and database\n";
+	return 1;
+    }
+
+
+    ////// INITIALISATION
+    // read Script and Database and put it in a std::string
+    auto script = readFile(argv[1]);
+    auto db     = readFile(argv[2]);
+
+    // Create the matrix directly from database and init script
+    std::cout << ":: Creating a PhreeqcMatrix with valence states and H(0)/O(0) \n";
+    PhreeqcMatrix pqc_mat1(db, script, true, true);
+    
+    // How many different SOLUTIONS ("CELLS") are defined in the script?
+    const auto ids = pqc_mat1.getIds();
+    
+    int n = ids.size();
+    
+    std::cout << ":: Found " << n << " distinct PHREEQC problems \n";
+
+    std::cout << ":: getSolutionsNames(): the common solutes across all problems: \n";
+    const auto solutes1  = pqc_mat1.getSolutionNames();
+    std::cout << solutes1 << "\n";
+
+    // auto expmat = pqc_mat1.get();
+    auto allvars = pqc_mat1.get().names;
+    std::cout << ":: pqc_mat1.get().names (all names in the PhreeqcMatrix): \n";
+    std::cout << allvars << "\n\n";
+
+    std::cout << "\n-- Now the new getMatrix*() --\n\n";
+    
+    auto transported = pqc_mat1.getMatrixTransported(); 
+    std::cout << ":: pqc_mat1.getMatrixTransported(): \n";
+    std::cout << transported << "\n\n";
+
+    auto MatNamesKin = pqc_mat1.getMatrixKinetics(); 
+    std::cout << ":: pqc_mat1.getMatrixKinetics(): \n";
+    std::cout << MatNamesKin << "\n\n";
+    
+    auto MatNamesEqui = pqc_mat1.getMatrixEquilibrium(); 
+    std::cout << ":: pqc_mat1.getMatrixEquilibrium(): \n";
+    std::cout << MatNamesEqui << "\n\n";
+
+    auto outonly = pqc_mat1.getMatrixOutOnly(); 
+    std::cout << ":: pqc_mat1.getMatrixOutOnly(): \n";
+    std::cout << outonly << "\n\n";
+    return 0;
+}
+
+// Oneliner for rz-vm278 relative to iphreeqc/poet/test!!
+
+// g++ testGolemRunner.cpp -o testG -Wall -I../../poet/include -I../../src -I../../src/phreeqcpp -I../../src/phreeqcpp/common -I../../src/phreeqcpp/PhreeqcKeywords -lIPhreeqc -lIPhreeqcPOET -L../../bbuild/ -L../../bbuild/poet

poet/test/testGolemRunner.cpp

@@ -1,4 +1,4 @@
-//  Time-stamp: "Last modified 2024-12-02 17:37:08 delucia"
+//  Time-stamp: "Last modified 2025-07-28 13:03:01 delucia"
 #include <iostream>
 #include <iomanip>
 #include <linux/limits.h>
@@ -60,7 +60,7 @@ int main(int argc, char *argv[]) {
     auto db     = readFile(argv[2]);
 
     // Create the matrix directly from database and init script
-    PhreeqcMatrix pqc_mat(db, script);
+    PhreeqcMatrix pqc_mat(db, script, true, true);
     
     // How many different SOLUTIONS ("CELLS") are defined in the script?
     const auto ids = pqc_mat.getIds();
@@ -99,12 +99,6 @@ int main(int argc, char *argv[]) {
     std::vector<double> &cell_values = exported_mat.values;
 
     std::cout << ":: Values in the PhreeqcMatrix: \n";
-
-    // std::cout << exported_mat.names  << "\n";
-    // std::cout << cell_values << "\n";
-    // END INIT
-
-
     
     //// Phreeqc RUN through the new Runner class 
 
@@ -117,39 +111,28 @@ int main(int argc, char *argv[]) {
     const auto matrix_values = stl_mat.values;
     const auto num_columns = stl_mat.names.size();
     const auto spec_names = stl_mat.names;
-
+    
     // container to pass in/out
     std::vector<std::vector<double>> simulationInOut;
 
-    // grid cells
-    const std::size_t num_cells = 10;
-    const std::size_t half_cells = 5;
-
     // copy the values to the InOut vector. We replicate cell 1 
-    for (std::size_t index = 0; index < num_cells; ++index) {
-	if (index < half_cells) {
-	    simulationInOut.push_back(std::vector<double>(
-            matrix_values.begin(), matrix_values.begin() + num_columns));
-	} else {
-	    simulationInOut.push_back(std::vector<double>(
-            matrix_values.begin() + num_columns, matrix_values.end()));
-	}
+    for (std::size_t index = 0; index < n; ++index) {
+	simulationInOut.push_back(std::vector<double>(
+        matrix_values.begin() + num_columns*index, matrix_values.begin() + num_columns*(index +1)));
     }
 
     const double timestep = 100.;
 
     // compute 1 timestep
     runner.run(simulationInOut, timestep);
-
+    
     for (std::size_t cell_index = 0; cell_index < simulationInOut.size(); ++cell_index) {
-	const bool is_first_half = cell_index < half_cells;
-	if (is_first_half) {
-	    std::cout << "Grid element: " << cell_index << " \n";
-	    for (std::size_t spec = 0; spec < num_columns; ++spec) {
-		std::cout << ":" << spec_names[spec] << "=" << simulationInOut[cell_index][spec];
-	    }
-	    std::cout << "\n";
+	
+	std::cout << "Grid element: " << cell_index << " \n";
+	for (std::size_t spec = 0; spec < num_columns; ++spec) {
+	    std::cout << ":" << spec_names[spec] << "=" << simulationInOut[cell_index][spec];
 	}
+	std::cout << "\n";
     }
 
     

src/phreeqcpp/GFZ.cpp

@@ -1,7 +1,8 @@
 #include "Phreeqc.h"
 #include <set>
 
-const std::set<std::string> to_ignore = {"H", "O", "Charge", "H(0)", "O(0)"};
+const std::set<std::string> to_ignore = {
+    "H", "O", "Charge", "tc", "patm", "SolVol", "pH", "pe", "H(0)", "O(0)"};
 
 std::vector<std::string> Phreeqc::find_all_valence_states(
     const std::vector<std::string> &solution_names) {
@@ -85,4 +86,4 @@ std::vector<std::string> Phreeqc::find_all_valence_states(
   }
 
   return solution_with_valences;
-}
+}

src/phreeqcpp/Solution.cxx

@@ -1061,26 +1061,26 @@ void cxxSolution::read_raw(CParser &parser, bool check) {
   return;
 }
 
-void cxxSolution::Update(LDBLE h_tot, LDBLE o_tot, LDBLE charge,
+void cxxSolution::Update(LDBLE h_tot, LDBLE o_tot, LDBLE charge, LDBLE tc, LDBLE patm,
                          const cxxNameDouble &const_nd) {
   this->new_def = false;
-  this->patm = 1.0;
-  this->potV = 0.0;
-  this->tc = 25.0;
-  this->ph = 7.0;
-  this->pe = 4.0;
-  this->mu = 1e-7;
-  this->ah2o = 1.0;
+  this->patm = patm;
+  // this->potV = 0.0;
+  this->tc = tc;
+  // this->ph = 7.0;
+  // this->pe = 4.0;
+  // this->mu = 1e-7;
+  // this->ah2o = 1.0;
   // H, O, charge, totals, and activities of solution are updated
   this->total_h = h_tot;
   this->total_o = o_tot;
   this->cb = charge;
   this->mass_water = o_tot / 55.55;
 
-  this->density = 1.0;
-  this->viscosity = 1.0;
-  this->soln_vol = 1.0;
-  this->total_alkalinity = 0.0;
+  // this->density = 1.0;
+  // this->viscosity = 1.0;
+  // this->soln_vol = 1.0;
+  // this->total_alkalinity = 0.0;
 
   this->master_activity.clear();
   this->species_gamma.clear();

src/phreeqcpp/Solution.h

@@ -123,7 +123,7 @@ public:
   // void modify_activities(const cxxSolution & original);
   // void Simplify_totals();
   void Update(const cxxNameDouble &nd);
-  void Update(LDBLE h_tot, LDBLE o_tot, LDBLE charge, const cxxNameDouble &nd);
+  void Update(LDBLE h_tot, LDBLE o_tot, LDBLE charge, LDBLE tc, LDBLE patm, const cxxNameDouble &nd);
   void Update_activities(const cxxNameDouble &original_tot);
   void Serialize(Dictionary &dictionary, std::vector<int> &ints,
                  std::vector<double> &doubles);