Engine: add error handling after run

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/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,9 +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_soln_vol();
-  data[4] = solution->Get_ph();
-  data[5] = solution->Get_pe();
+  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()
@@ -44,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++];
@@ -54,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());
 }

poet/src/Wrapper/SolutionWrapper.hpp

@@ -27,9 +27,10 @@ private:
   cxxSolution *solution;
   const std::vector<std::string> solution_order;
 
-  static constexpr std::array<const char *, 6> ESSENTIALS = {"H", "O",
-                                                             "Charge",
-                                                             "SolVol", "pH","pe"}; // MDL
+  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();
 

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/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);