BREAKING CHANGE: Add wrappers for reactants to get/set values

poet/include/EquilibriumWrapper.hpp

@@ -0,0 +1,38 @@
+#pragma once
+
+#include "PPassemblage.h"
+#include "WrapperBase.hpp"
+#include <memory>
+
+class EquilibriumWrapper : public WrapperBase {
+public:
+  EquilibriumWrapper(cxxPPassemblage *ppassemblage,
+                     const std::vector<std::string> &ppassemblage_names);
+
+  void get(std::span<LDBLE> &data) const override;
+
+  void set(const std::span<LDBLE> &data) override;
+
+  static std::vector<std::string>
+  names(const cxxPPassemblage *ppassemblage,
+        std::vector<std::string> &ppassemblage_names);
+
+private:
+  cxxPPassemblage *ppassemblage;
+
+  class EquilibriumCompWrapper : public WrapperBase {
+  public:
+    EquilibriumCompWrapper(cxxPPassemblageComp &comp);
+
+    void get(std::span<LDBLE> &data) const override;
+
+    void set(const std::span<LDBLE> &data) override;
+
+    static std::vector<std::string> names(const cxxPPassemblageComp &comp);
+
+  private:
+    cxxPPassemblageComp &comp;
+  };
+
+  std::vector<std::unique_ptr<EquilibriumCompWrapper>> equilibrium_comps;
+};

poet/include/ExchangeWrapper.hpp

@@ -0,0 +1,43 @@
+#pragma once
+
+#include "ExchComp.h"
+#include "Exchange.h"
+#include "WrapperBase.hpp"
+#include "phrqtype.h"
+#include <cstddef>
+#include <memory>
+#include <string>
+#include <vector>
+
+#include <span>
+
+class ExchangeWrapper : public WrapperBase {
+public:
+  ExchangeWrapper(cxxExchange *exch, const std::vector<std::string> &exchanger);
+
+  void get(std::span<LDBLE> &exchange) const;
+
+  void set(const std::span<LDBLE> &exchange);
+
+  static std::vector<std::string>
+  names(cxxExchange *exchange, std::vector<std::string> &exchange_formulas);
+
+private:
+  cxxExchange *exchange;
+
+  class ExchangeCompWrapper : public WrapperBase {
+  public:
+    ExchangeCompWrapper(cxxExchComp &comp);
+
+    void get(std::span<LDBLE> &exchange) const;
+
+    void set(const std::span<LDBLE> &exchange);
+
+  private:
+    cxxExchComp &exch_comp;
+
+    static constexpr std::size_t NUM_NOT_TOTALS = 5;
+  };
+
+  std::vector<std::unique_ptr<ExchangeCompWrapper>> exchange_comps;
+};

poet/include/IPhreeqcPOET.hpp

@@ -1,192 +0,0 @@
-#pragma once
-
-#include <IPhreeqc.hpp>
-
-#include <Exchange.h>
-#include <PPassemblage.h>
-#include <Phreeqc.h>
-#include <Solution.h>
-#include <Surface.h>
-#include <array>
-#include <cstddef>
-#include <cstdint>
-#include <cxxKinetics.h>
-#include <iomanip>
-#include <map>
-#include <sstream>
-#include <string>
-#include <sys/types.h>
-#include <utility>
-#include <vector>
-
-enum { POET_SOL = 0, POET_EXCH, POET_KIN, POET_EQUIL, POET_SURF };
-
-class IPhreeqcPOET : public IPhreeqc {
-public:
-  IPhreeqcPOET(const std::string &database, const std::string &input_script)
-      : IPhreeqc() {
-    this->LoadDatabaseString(database.c_str());
-    this->RunString(input_script.c_str());
-
-    this->parseInit();
-  }
-
-  IPhreeqcPOET(const std::string &database, const std::string &input_script,
-               const std::vector<std::string> &solutionInitVector,
-               std::uint32_t n_cells)
-      : IPhreeqc(), n_cells(n_cells), solutionInitVector(solutionInitVector) {
-    this->LoadDatabaseString(database.c_str());
-    this->RunString(input_script.c_str());
-
-    if (n_cells > 1) {
-      std::string copy_string =
-          "COPY cell 1 2-" + std::to_string(n_cells) + "\n";
-      this->RunString(copy_string.c_str());
-    }
-  }
-
-  void queueCell(std::vector<double> &values) {
-    if (this->queued_cell_pointer > this->n_cells) {
-      return;
-    }
-
-    setCell(this->queued_cell_pointer++, values);
-  }
-
-  void dequeueCells(std::vector<std::vector<double>> &values) {
-    if (this->queued_cell_pointer == 1) {
-      return;
-    }
-
-    for (std::size_t i = 1; i < this->queued_cell_pointer; i++) {
-      std::vector<double> cell_values;
-      getCell(i, cell_values);
-      values.push_back(cell_values);
-    }
-
-    this->queued_cell_pointer = 1;
-  }
-
-  void runQueuedCells(double time_step) {
-    if (this->queued_cell_pointer == 1) {
-      return;
-    }
-    run(1, this->queued_cell_pointer - 1, time_step);
-  }
-
-  void setCell(int cell_number, std::vector<double> &values) {
-    this->set_essential_values(cell_number, this->solutionInitVector, values);
-  }
-
-  void getCell(int cell_number, std::vector<double> &values) {
-    values = this->get_essential_values(cell_number, this->solutionInitVector);
-  }
-
-  void run(int start_cell, int end_cell, double time_step) {
-    std::stringstream time_ss;
-    time_ss << std::fixed << std::setprecision(20) << time_step;
-
-    const std::string runs_string =
-        "RUN_CELLS\n -cells " + std::to_string(start_cell) + "-" +
-        std::to_string(end_cell) + "\n -time_step " + time_ss.str();
-    this->RunString(runs_string.c_str());
-  }
-
-  struct PhreeqcMat {
-    std::vector<std::string> names;
-    std::vector<int> ids;
-    std::vector<std::vector<double>> values;
-  };
-
-  PhreeqcMat getPhreeqcMat();
-
-  std::map<int, std::string> raw_dumps() {
-    std::map<int, std::string> 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);
-      this->RunString(call_string.c_str());
-      dumps[sol_id] = this->GetDumpString();
-    }
-
-    this->SetDumpStringOn(false);
-
-    return dumps;
-  }
-
-  using essential_names = std::array<std::vector<std::string>, 5>;
-  using ModulesArray = std::array<std::uint32_t, 5>;
-
-  ModulesArray getModuleSizes(const std::vector<int> &cell_ids);
-
-  std::vector<std::string> getSolutionNames(int cell_id) {
-    return this->raw_initials[cell_id].first[POET_SOL];
-  }
-
-private:
-  // required only for simulation
-  essential_names initial_names;
-  std::uint32_t n_cells;
-  std::vector<std::string> solutionInitVector;
-
-  void parseInitValues();
-  void parseInit();
-
-  essential_names dump_essential_names(std::size_t cell_number);
-
-  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);
-  }
-
-  using RawMap = std::map<int, std::pair<essential_names, std::vector<double>>>;
-
-  essential_names union_raws(const RawMap &raws);
-
-  std::vector<std::vector<double>>
-  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<double> &values);
-
-  std::string 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;
-  };
-
-  std::vector<double>
-  get_essential_values(std::size_t cell_number,
-                       const std::vector<std::string> &order);
-
-  void set_essential_values(std::size_t cell_number,
-                            const std::vector<std::string> &order,
-                            std::vector<double> &values);
-
-  RawMap raw_initials;
-
-  std::size_t queued_cell_pointer = 1;
-};

poet/include/KineticWrapper.hpp

@@ -0,0 +1,38 @@
+#pragma once
+
+#include "KineticsComp.h"
+#include "WrapperBase.hpp"
+#include "cxxKinetics.h"
+#include <memory>
+#include <vector>
+
+class KineticWrapper : public WrapperBase {
+
+public:
+  KineticWrapper(cxxKinetics *kinetics,
+                 const std::vector<std::string> &kin_comps);
+
+  void get(std::span<LDBLE> &data) const override;
+  void set(const std::span<LDBLE> &data) override;
+
+  static std::vector<std::string> names(cxxKinetics *kinetics,
+                                        std::vector<std::string> &kin_comps);
+
+private:
+  cxxKinetics *kinetics;
+
+  class KineticCompWrapper : public WrapperBase {
+  public:
+    KineticCompWrapper(cxxKineticsComp &comp);
+
+    void get(std::span<LDBLE> &data) const override;
+    void set(const std::span<LDBLE> &data) override;
+
+    static std::vector<std::string> names(const cxxKineticsComp &comp);
+
+  private:
+    cxxKineticsComp &kin_comp;
+  };
+
+  std::vector<std::unique_ptr<KineticCompWrapper>> kinetic_comps;
+};

poet/include/POETInit.hpp

@@ -0,0 +1,20 @@
+#pragma once
+
+#include <string>
+#include <vector>
+
+struct POETInitCell {
+  std::vector<std::string> solutions;
+  std::vector<std::string> solution_primaries;
+  std::vector<std::string> exchanger;
+  std::vector<std::string> kinetics;
+  std::vector<std::string> equilibrium;
+  std::vector<std::string> surface_comps;
+  std::vector<std::string> surface_charges;
+};
+
+struct POETConfig {
+  std::string database;
+  std::string input_script;
+  POETInitCell cell;
+};

poet/include/PhreeqcEngine.hpp

@@ -0,0 +1,96 @@
+
+#pragma once
+
+#include "EquilibriumWrapper.hpp"
+#include "ExchangeWrapper.hpp"
+#include "KineticWrapper.hpp"
+#include "POETInit.hpp"
+#include "SolutionWrapper.hpp"
+#include "SurfaceWrapper.hpp"
+
+#include <IPhreeqc.hpp>
+
+#include <Exchange.h>
+#include <PPassemblage.h>
+#include <Phreeqc.h>
+#include <Solution.h>
+#include <Surface.h>
+#include <cstddef>
+#include <cxxKinetics.h>
+#include <iomanip>
+#include <iostream>
+#include <memory>
+#include <ostream>
+#include <span>
+#include <sstream>
+#include <string>
+#include <sys/types.h>
+#include <vector>
+
+class PhreeqcEngine : public IPhreeqc {
+public:
+  PhreeqcEngine(const POETConfig &config) : IPhreeqc() {
+    this->LoadDatabaseString(config.database.c_str());
+    this->RunString(config.input_script.c_str());
+
+    this->init_wrappers(config.cell);
+  }
+
+  void runCell(std::vector<double> &cell_values, double time_step);
+
+  void run(int start_cell, int end_cell) {
+    const std::string runs_string = "RUN_CELLS\n -cells " +
+                                    std::to_string(start_cell) + "-" +
+                                    std::to_string(end_cell) + "\nEND\n";
+    this->RunString(runs_string.c_str());
+  }
+
+  void run(int start_cell, int end_cell, double time_step) {
+    std::stringstream time_ss;
+    time_ss << std::fixed << std::setprecision(20) << time_step;
+
+    const std::string runs_string =
+        "RUN_CELLS\n -cells " + std::to_string(start_cell) + "-" +
+        std::to_string(end_cell) + "\n -time_step " + time_ss.str() + "\nEND\n";
+    this->RunString(runs_string.c_str());
+  }
+
+private:
+  void init_wrappers(const POETInitCell &cell);
+
+  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);
+  }
+
+  void get_essential_values(std::span<double> &data);
+
+  void set_essential_values(const std::span<double> &data);
+
+  std::unique_ptr<SolutionWrapper> solutionWrapperPtr;
+  std::unique_ptr<ExchangeWrapper> exchangeWrapperPtr;
+  std::unique_ptr<KineticWrapper> kineticsWrapperPtr;
+  std::unique_ptr<EquilibriumWrapper> equilibriumWrapperPtr;
+  std::unique_ptr<SurfaceWrapper> surfaceWrapperPtr;
+
+  bool has_exchange = false;
+  bool has_kinetics = false;
+  bool has_equilibrium = false;
+  bool has_surface = false;
+};

poet/include/PhreeqcInit.hpp

@@ -0,0 +1,156 @@
+#pragma once
+
+#include <IPhreeqc.hpp>
+
+#include <Phreeqc.h>
+#include <array>
+#include <cstddef>
+#include <cstdint>
+#include <map>
+#include <string>
+#include <utility>
+#include <vector>
+
+#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 };
+
+class PhreeqcInit : public IPhreeqc {
+public:
+  PhreeqcInit(const std::string &database, const std::string &input_script);
+
+  struct PhreeqcMat {
+    std::vector<std::string> names;
+    std::vector<int> ids;
+    std::vector<std::vector<double>> values;
+  };
+
+  PhreeqcMat getPhreeqcMat() const { return pqc_mat; };
+
+  std::map<int, std::string> raw_dumps();
+
+  using essential_names = std::array<std::vector<std::string>, 5>;
+  using ModulesArray = std::array<std::uint32_t, 5>;
+
+  ModulesArray getModuleSizes(const std::vector<int> &cell_ids);
+
+  std::vector<std::string> getSolutionPrimaries() {
+    return std::vector<std::string>(this->surface_primaries.begin(),
+                                    this->surface_primaries.end());
+  }
+
+  std::vector<std::string> getSolutionNames(int cell_id) {
+    return this->raw_initials[cell_id].first[POET_SOL];
+  }
+
+  std::vector<std::string> getExchanger(int id) {
+    if (this->exchanger.contains(id)) {
+      return this->exchanger[id];
+    }
+
+    return {};
+  }
+
+  std::vector<std::string> getKineticsNames(int id) {
+    if (this->kinetics.contains(id)) {
+      return this->kinetics[id];
+    }
+
+    return {};
+  }
+
+  std::vector<std::string> getEquilibriumNames(int id) {
+    if (this->equilibrium.contains(id)) {
+      return this->equilibrium[id];
+    }
+
+    return {};
+  }
+
+  std::vector<std::string> getSurfaceCompNames(int id) {
+    if (this->surface_comps.contains(id)) {
+      return this->surface_comps[id];
+    }
+
+    return {};
+  }
+
+  std::vector<std::string> getSurfaceChargeNames(int id) {
+    if (this->surface_charge.contains(id)) {
+      return this->surface_charge[id];
+    }
+
+    return {};
+  }
+
+private:
+  // required only for simulation
+  essential_names initial_names;
+
+  void parseInitValues();
+  void parseInit();
+
+  std::vector<std::string> dump_solution_names(int cell_number);
+  void dump_reactant_names(int cell_number,
+                           const std::vector<std::string> union_sol_names,
+                           essential_names &names);
+
+  std::vector<std::string>
+  find_all_valence_states(const std::vector<std::string> &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<int, std::vector<std::string>> exchanger;
+  std::map<int, std::vector<std::string>> kinetics;
+  std::map<int, std::vector<std::string>> equilibrium;
+  std::map<int, std::vector<std::string>> surface_comps;
+  std::map<int, std::vector<std::string>> surface_charge;
+
+  std::set<std::string> surface_primaries;
+
+  using RawMap = std::map<int, std::pair<essential_names, std::vector<double>>>;
+
+  essential_names union_raws(const RawMap &raws);
+
+  std::vector<std::vector<double>>
+  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<double> &values);
+
+  std::string subExchangeName(std::string name);
+
+  std::vector<double>
+  get_essential_values_init(std::size_t cell_number,
+                            const std::vector<std::string> &order);
+
+  RawMap raw_initials;
+};

poet/include/SolutionWrapper.hpp

@@ -0,0 +1,26 @@
+#pragma once
+
+#include "NameDouble.h"
+#include "Solution.h"
+#include "WrapperBase.hpp"
+#include <cstddef>
+#include <vector>
+
+class SolutionWrapper : public WrapperBase {
+public:
+  SolutionWrapper(cxxSolution *soln,
+                  const std::vector<std::string> &solution_order);
+
+  void get(std::span<LDBLE> &data) const;
+
+  void set(const std::span<LDBLE> &data);
+
+  static std::vector<std::string>
+  names(cxxSolution *solution, std::vector<std::string> &solution_order);
+
+private:
+  cxxSolution *solution;
+  const std::vector<std::string> solution_order;
+
+  static constexpr std::size_t NUM_ESSENTIALS = 3;
+};

poet/include/SurfaceWrapper.hpp

@@ -0,0 +1,71 @@
+#pragma once
+
+#include "Surface.h"
+#include "SurfaceCharge.h"
+#include "SurfaceComp.h"
+#include "WrapperBase.hpp"
+#include "phrqtype.h"
+#include <cstddef>
+#include <memory>
+#include <set>
+#include <span>
+#include <string>
+#include <vector>
+
+class SurfaceWrapper : public WrapperBase {
+public:
+  SurfaceWrapper(cxxSurface *surf,
+                 const std::set<std::string> &solution_primaries,
+                 const std::vector<std::string> &comp_formulas,
+                 const std::vector<std::string> &charge_names);
+
+  void get(std::span<LDBLE> &surface) const override;
+
+  void set(const std::span<LDBLE> &surface) override;
+
+  static std::vector<std::string>
+  names(cxxSurface *surface, const std::set<std::string> &solution_primaries,
+        std::vector<std::string> &comp_formulas,
+        std::vector<std::string> &charge_names);
+
+private:
+  cxxSurface *surface;
+
+  class SurfaceCompWrapper : public WrapperBase {
+  public:
+    SurfaceCompWrapper(cxxSurfaceComp &comp);
+
+    void get(std::span<LDBLE> &surface) const;
+
+    void set(const std::span<LDBLE> &surface);
+
+    static std::vector<std::string> names(cxxSurfaceComp &comp);
+
+  private:
+    cxxSurfaceComp &surface_comp;
+    static constexpr std::size_t NUM_NOT_TOTALS = 3;
+
+    std::vector<std::string> total_names;
+  };
+
+  class SurfaceChargeWrapper : public WrapperBase {
+  public:
+    SurfaceChargeWrapper(cxxSurfaceCharge &charge,
+                         const std::set<std::string> &solution_primaries);
+    void get(std::span<LDBLE> &surface) const;
+
+    void set(const std::span<LDBLE> &surface);
+
+    static std::vector<std::string>
+    names(cxxSurfaceCharge *charge, const std::set<std::string> &primaries);
+
+  private:
+    cxxSurfaceCharge &surface_charge;
+    static constexpr std::size_t NUM_NOT_TOTALS = 5;
+
+    std::set<std::string> primaries;
+  };
+
+  std::vector<std::unique_ptr<SurfaceCompWrapper>> surface_comps;
+  std::vector<std::unique_ptr<SurfaceChargeWrapper>> surface_charges;
+};

poet/include/WrapperBase.hpp

@@ -0,0 +1,18 @@
+#pragma once
+
+#include <phrqtype.h>
+#include <span>
+
+class WrapperBase {
+public:
+  virtual ~WrapperBase() = default;
+
+  std::size_t size() const { return this->num_elements; };
+
+  virtual void get(std::span<LDBLE> &data) const = 0;
+
+  virtual void set(const std::span<LDBLE> &data) = 0;
+
+protected:
+  std::size_t num_elements = 0;
+};

poet/src/Engine.cpp

@@ -0,0 +1,220 @@
+#include "PhreeqcEngine.hpp"
+#include <cstddef>
+#include <span>
+
+void PhreeqcEngine::init_wrappers(const POETInitCell &cell) {
+
+  // TODO: Implement the rest of the wrappers. Currently supported: EXCHANGE
+
+  // Solutions
+  this->solutionWrapperPtr =
+      std::make_unique<SolutionWrapper>(this->Get_solution(1), cell.solutions);
+
+  if (this->Get_exchange(1) != nullptr) {
+    this->exchangeWrapperPtr = std::make_unique<ExchangeWrapper>(
+        this->Get_exchange(1), cell.exchanger);
+    this->has_exchange = true;
+  }
+
+  if (this->Get_kinetic(1) != nullptr) {
+    this->kineticsWrapperPtr =
+        std::make_unique<KineticWrapper>(this->Get_kinetic(1), cell.kinetics);
+
+    this->has_kinetics = true;
+  }
+
+  if (this->Get_equilibrium(1) != nullptr) {
+    this->equilibriumWrapperPtr = std::make_unique<EquilibriumWrapper>(
+        this->Get_equilibrium(1), cell.equilibrium);
+
+    this->has_equilibrium = true;
+  }
+
+  if (this->Get_surface(1) != nullptr) {
+    std::set<std::string> primaries(cell.solution_primaries.begin(),
+                                    cell.solution_primaries.end());
+    this->surfaceWrapperPtr = std::make_unique<SurfaceWrapper>(
+        this->Get_surface(1), primaries, cell.surface_comps,
+        cell.surface_charges);
+
+    this->has_surface = true;
+  }
+}
+
+void PhreeqcEngine::get_essential_values(std::span<double> &data) {
+
+  this->solutionWrapperPtr->get(data);
+
+  std::size_t offset = this->solutionWrapperPtr->size();
+
+  if (this->has_exchange) {
+    std::span<double> exch_span{
+        data.subspan(offset, this->exchangeWrapperPtr->size())};
+    this->exchangeWrapperPtr->get(exch_span);
+
+    offset += this->exchangeWrapperPtr->size();
+  }
+
+  if (this->has_kinetics) {
+    std::span<double> kin_span{
+        data.subspan(offset, this->kineticsWrapperPtr->size())};
+    this->kineticsWrapperPtr->get(kin_span);
+
+    offset += this->kineticsWrapperPtr->size();
+  }
+
+  if (this->has_equilibrium) {
+    std::span<double> equ_span{
+        data.subspan(offset, this->equilibriumWrapperPtr->size())};
+    this->equilibriumWrapperPtr->get(equ_span);
+
+    offset += this->equilibriumWrapperPtr->size();
+  }
+
+  if (this->has_surface) {
+    std::span<double> surf_span{
+        data.subspan(offset, this->surfaceWrapperPtr->size())};
+    this->surfaceWrapperPtr->get(surf_span);
+  }
+
+  // // Exchange
+  // if (this->Get_exchange(cell_number) != NULL) {
+  //   std::vector<double> exch_values(this->exchangeWrapperPtr->size());
+
+  //   std::span<double> exch_span{exch_values};
+
+  //   this->exchangeWrapperPtr->get(exch_span);
+
+  //   // this->Get_exchange(cell_number)->get_essential_values(exch_values);
+  //   essentials.insert(essentials.end(), exch_values.begin(),
+  //   exch_values.end());
+  // }
+
+  // // Kinetics
+  // if (this->Get_kinetic(cell_number) != NULL) {
+  //   std::vector<double> kin_values;
+
+  //   this->Get_kinetic(cell_number)->get_essential_values(kin_values);
+  //   essentials.insert(essentials.end(), kin_values.begin(),
+  //   kin_values.end());
+  // }
+
+  // // PPassemblage
+  // if (this->Get_equilibrium(cell_number) != NULL) {
+  //   std::vector<double> equ_values;
+
+  //   this->Get_equilibrium(cell_number)->get_essential_values(equ_values);
+  //   essentials.insert(essentials.end(), equ_values.begin(),
+  //   equ_values.end());
+  // }
+
+  // // Surface
+  // if (this->Get_surface(cell_number) != NULL) {
+  //   std::vector<double> surf_values(this->surfaceWrapperPtr->size());
+  //   std::span<double> surf_span{surf_values};
+
+  //   this->surfaceWrapperPtr->get(surf_span);
+
+  //   // this->Get_surface(cell_number)->get_essential_values(surf_values);
+  //   essentials.insert(essentials.end(), surf_values.begin(),
+  //   surf_values.end());
+  // }
+
+  // return essentials;
+}
+
+void PhreeqcEngine::set_essential_values(const std::span<double> &data) {
+
+  this->solutionWrapperPtr->set(data);
+  this->PhreeqcPtr->initial_solutions_poet(1);
+
+  std::size_t offset = this->solutionWrapperPtr->size();
+
+  if (this->has_exchange) {
+    std::span<double> exch_span{
+        data.subspan(offset, this->exchangeWrapperPtr->size())};
+    this->exchangeWrapperPtr->set(exch_span);
+
+    offset += this->exchangeWrapperPtr->size();
+  }
+
+  if (this->has_kinetics) {
+    std::span<double> kin_span{
+        data.subspan(offset, this->kineticsWrapperPtr->size())};
+    this->kineticsWrapperPtr->set(kin_span);
+
+    offset += this->kineticsWrapperPtr->size();
+  }
+
+  if (this->has_equilibrium) {
+    std::span<double> equ_span{
+        data.subspan(offset, this->equilibriumWrapperPtr->size())};
+    this->equilibriumWrapperPtr->set(equ_span);
+
+    offset += this->equilibriumWrapperPtr->size();
+  }
+
+  if (this->has_surface) {
+    std::span<double> surf_span{
+        data.subspan(offset, this->surfaceWrapperPtr->size())};
+    this->surfaceWrapperPtr->set(surf_span);
+  }
+
+  // // Solutions
+  // if (this->Get_solution(cell_number) != NULL) {
+  //   this->Get_solution(cell_number)->set_essential_values(dump_it, order);
+  //   this->PhreeqcPtr->initial_solutions_poet(cell_number);
+  // }
+
+  // // Exchange
+  // if (this->Get_exchange(cell_number) != NULL) {
+  //   auto &exch_pointer = this->exchangeWrapperPtr;
+  //   std::span<double> exch_span{dump_it, exch_pointer->size()};
+  //   exch_pointer->set(exch_span);
+
+  //   dump_it += exch_pointer->size();
+
+  //   // this->Get_exchange(cell_number)->set_essential_values(dump_it);
+  //   // this->PhreeqcPtr->Rxn_new_exchange.insert(cell_number);
+  //   // this->Get_exchange(cell_number)->Set_new_def(true);
+  //   // this->Get_exchange(cell_number)->Set_solution_equilibria(true);
+  //   // this->Get_exchange(cell_number)->Set_n_solution(cell_number);
+  //   // this->PhreeqcPtr->initial_exchangers(FALSE);
+  // }
+
+  // // Kinetics
+  // if (this->Get_kinetic(cell_number) != NULL) {
+  //   this->Get_kinetic(cell_number)->set_essential_values(dump_it);
+  // }
+
+  // // PPassemblage
+  // if (this->Get_equilibrium(cell_number) != NULL) {
+  //   this->Get_equilibrium(cell_number)->set_essential_values(dump_it);
+  // }
+
+  // // Surface
+  // if (this->Get_surface(cell_number) != NULL) {
+  //   auto *tmp = this->Get_surface(cell_number);
+  //   auto &surf_pointer = this->surfaceWrapperPtr;
+  //   std::span<double> surf_span{dump_it, surf_pointer->size()};
+  //   surf_pointer->set(surf_span);
+
+  //   dump_it += surf_pointer->size();
+
+  //   // this->PhreeqcPtr->Rxn_new_surface.insert(cell_number);
+
+  //   // this->Get_surface(cell_number)->set_essential_values(dump_it);
+  //   // this->PhreeqcPtr->Rxn_new_surface.insert(cell_number);
+  //   // this->PhreeqcPtr->initial_surfaces(TRUE);
+  // }
+}
+
+void PhreeqcEngine::runCell(std::vector<double> &cell_values,
+                            double time_step) {
+  // skip ID
+  std::span<double> cell_data{cell_values.begin() + 1, cell_values.end()};
+
+  this->set_essential_values(cell_data);
+  this->run(1, 1, time_step);
+  this->get_essential_values(cell_data);
+}

poet/src/EquilibriumCompWrapper.cpp

@@ -0,0 +1,30 @@
+#include "EquilibriumWrapper.hpp"
+
+EquilibriumWrapper::EquilibriumCompWrapper::EquilibriumCompWrapper(
+    cxxPPassemblageComp &comp_)
+    : comp(comp_) {
+  this->num_elements = 2;
+}
+
+void EquilibriumWrapper::EquilibriumCompWrapper::get(
+    std::span<LDBLE> &data) const {
+  data[0] = this->comp.Get_moles();
+  data[1] = this->comp.Get_si();
+}
+
+void EquilibriumWrapper::EquilibriumCompWrapper::set(
+    const std::span<LDBLE> &data) {
+  this->comp.Set_moles(data[0]);
+  this->comp.Set_si(data[1]);
+}
+
+std::vector<std::string> EquilibriumWrapper::EquilibriumCompWrapper::names(
+    const cxxPPassemblageComp &comp) {
+  std::vector<std::string> names;
+
+  const std::string &comp_name = comp.Get_name();
+  names.push_back(comp_name);
+  names.push_back(comp_name + "_si");
+
+  return names;
+}

poet/src/EquilibriumWrapper.cpp

@@ -0,0 +1,66 @@
+#include "EquilibriumWrapper.hpp"
+
+#include <algorithm>
+#include <memory>
+
+EquilibriumWrapper::EquilibriumWrapper(
+    cxxPPassemblage *ppassemblage,
+    const std::vector<std::string> &ppassemblage_names)
+    : ppassemblage(ppassemblage) {
+  for (const auto &ppassemblage_name : ppassemblage_names) {
+    auto it = std::find_if(
+        ppassemblage->Get_pp_assemblage_comps().begin(),
+        ppassemblage->Get_pp_assemblage_comps().end(),
+        [&](const auto &comp) { return comp.first == ppassemblage_name; });
+
+    if (it == ppassemblage->Get_pp_assemblage_comps().end()) {
+      throw std::runtime_error(
+          "Equilibrium component not found in Phreeqc variables");
+    }
+
+    equilibrium_comps.push_back(
+        std::make_unique<EquilibriumCompWrapper>(it->second));
+
+    num_elements += equilibrium_comps.back()->size();
+  }
+}
+
+void EquilibriumWrapper::get(std::span<LDBLE> &data) const {
+  std::size_t offset = 0;
+
+  for (const auto &comp : equilibrium_comps) {
+    std::span<LDBLE> comp_span = data.subspan(offset, comp->size());
+
+    comp->get(comp_span);
+
+    offset += comp->size();
+  }
+}
+
+void EquilibriumWrapper::set(const std::span<LDBLE> &data) {
+  std::size_t offset = 0;
+
+  for (const auto &comp : equilibrium_comps) {
+    std::span<LDBLE> comp_span = data.subspan(offset, comp->size());
+
+    comp->set(comp_span);
+
+    offset += comp->size();
+  }
+}
+
+std::vector<std::string>
+EquilibriumWrapper::names(const cxxPPassemblage *ppassemblage,
+                          std::vector<std::string> &ppassemblage_names) {
+  std::vector<std::string> names;
+
+  for (const auto &comp : ppassemblage->Get_pp_assemblage_comps()) {
+    std::vector<std::string> comp_names =
+        EquilibriumCompWrapper::names(comp.second);
+    names.insert(names.end(), comp_names.begin(), comp_names.end());
+
+    ppassemblage_names.push_back(comp.first);
+  }
+
+  return names;
+}

poet/src/ExchangeCompWrapper.cpp

@@ -0,0 +1,41 @@
+#include "ExchangeWrapper.hpp"
+
+ExchangeWrapper::ExchangeCompWrapper::ExchangeCompWrapper(cxxExchComp &comp)
+    : exch_comp(comp) {
+  const std::size_t totals_size = exch_comp.Get_totals().size() - 1;
+  this->num_elements = totals_size + NUM_NOT_TOTALS;
+}
+
+void ExchangeWrapper::ExchangeCompWrapper::get(
+    std::span<LDBLE> &exchange) const {
+  exchange[0] = exch_comp.Get_totals().find(exch_comp.Get_formula())->second;
+  exchange[1] = exch_comp.Get_charge_balance();
+  exchange[2] = exch_comp.Get_la();
+  exchange[3] = exch_comp.Get_phase_proportion();
+  exchange[4] = exch_comp.Get_formula_z();
+
+  for (const auto &[name, value] : exch_comp.Get_totals()) {
+    if (name != exch_comp.Get_formula()) {
+      exchange[NUM_NOT_TOTALS +
+               std::distance(exch_comp.Get_totals().begin(),
+                             exch_comp.Get_totals().find(name))] = value;
+    }
+  }
+}
+
+void ExchangeWrapper::ExchangeCompWrapper::set(
+    const std::span<LDBLE> &exchange) {
+  exch_comp.Get_totals().find(exch_comp.Get_formula())->second = exchange[0];
+  exch_comp.Set_charge_balance(exchange[1]);
+  exch_comp.Set_la(exchange[2]);
+  exch_comp.Set_phase_proportion(exchange[3]);
+  exch_comp.Set_formula_z(exchange[4]);
+
+  for (auto &[name, value] : exch_comp.Get_totals()) {
+    if (name != exch_comp.Get_formula()) {
+      value = exchange[NUM_NOT_TOTALS +
+                       std::distance(exch_comp.Get_totals().begin(),
+                                     exch_comp.Get_totals().find(name))];
+    }
+  }
+}

poet/src/ExchangeWrapper.cpp

@@ -0,0 +1,101 @@
+#include "ExchangeWrapper.hpp"
+
+ExchangeWrapper::ExchangeWrapper(cxxExchange *exch,
+                                 const std::vector<std::string> &exchanger)
+    : exchange(exch) {
+
+  for (const auto &comp_name : exchanger) {
+    auto it = std::find_if(exchange->Get_exchange_comps().begin(),
+                           exchange->Get_exchange_comps().end(),
+                           [&](const cxxExchComp &comp) {
+                             return comp.Get_formula() == comp_name;
+                           });
+
+    if (it == exchange->Get_exchange_comps().end()) {
+      throw std::runtime_error(
+          "Exchange component not found in Phreeqc variables");
+    }
+
+    exchange_comps.push_back(std::make_unique<ExchangeCompWrapper>(*it));
+
+    num_elements += exchange_comps.back()->size();
+  }
+
+    // const std::size_t defined_comps = exchange->Get_exchange_comps().size();
+
+    // auto header_it = remaining_field_header.begin();
+
+    // while (header_it != remaining_field_header.end() &&
+    //        exchange_comps.size() < defined_comps) {
+    //   const std::string formular = *header_it;
+
+    //   auto it = std::find_if(exchange->Get_exchange_comps().begin(),
+    //                          exchange->Get_exchange_comps().end(),
+    //                          [&](const cxxExchComp &comp) {
+    //                            return comp.Get_formula() == formular;
+    //                          });
+
+    //   if (it != exchange->Get_exchange_comps().end()) {
+    //     const size_t i = this->exchange_comps.size();
+
+    //     exchange_comps.push_back(std::make_unique<ExchangeCompWrapper>(*it));
+    //     header_it += this->exchange_comps[i]->size();
+    //     num_elements += this->exchange_comps[i]->size();
+    //     continue;
+    //   }
+
+    //   header_it++;
+    // }
+
+    // if (exchange_comps.size() != defined_comps) {
+    //   throw std::runtime_error(
+    //       "Not all exchange components found in Phreeqc variables");
+    // }
+}
+
+void ExchangeWrapper::get(std::span<LDBLE> &exchange) const {
+  std::size_t offset = 0;
+
+  for (const auto &comp : exchange_comps) {
+    std::span<LDBLE> comp_span = exchange.subspan(offset, comp->size());
+    comp->get(comp_span);
+    offset += comp->size();
+  }
+}
+
+void ExchangeWrapper::set(const std::span<LDBLE> &exchange) {
+  std::size_t offset = 0;
+
+  for (const auto &comp : exchange_comps) {
+    std::span<LDBLE> comp_span = exchange.subspan(offset, comp->size());
+    comp->set(comp_span);
+    offset += comp->size();
+  }
+}
+
+std::vector<std::string>
+ExchangeWrapper::names(cxxExchange *exchange,
+                       std::vector<std::string> &exchange_formulas) {
+  std::vector<std::string> e_names;
+
+  for (const auto &comp : exchange->Get_exchange_comps()) {
+
+    const std::string &formular = comp.Get_formula();
+    exchange_formulas.push_back(formular);
+    e_names.push_back(formular);
+
+    e_names.push_back(formular + "_cb");
+    e_names.push_back(formular + "_la");
+    e_names.push_back(formular + "_phase_proportion");
+    e_names.push_back(formular + "_formular_z");
+
+    for (const auto &total : comp.Get_totals()) {
+      if (total.first == formular) {
+        continue;
+      }
+      e_names.push_back(total.first + formular);
+    }
+  }
+
+  return e_names;
+}

poet/src/Extraction.cpp

@@ -1,4 +1,6 @@
-#include <IPhreeqcPOET.hpp>
+#include <PhreeqcInit.hpp>
+
+#include "Solution.h"
 
 #include <algorithm>
 #include <cstddef>
@@ -37,9 +39,19 @@ createConcVector(const std::vector<std::string> &conc_names,
   return conc_vec;
 }
 
-void IPhreeqcPOET::parseInit() {
+PhreeqcInit::PhreeqcInit(const std::string &database,
+                         const std::string &input_script)
+    : IPhreeqc() {
+  this->LoadDatabaseString(database.c_str());
+  this->RunString(input_script.c_str());
 
-  std::map<int, std::pair<essential_names, std::vector<double>>> init_values;
+  this->parseInit();
+  pqc_mat = this->buildPhreeqcMat();
+}
+
+void PhreeqcInit::parseInit() {
+
+  std::map<int, essential_names> init_values;
 
   for (const auto &[id, val] : this->PhreeqcPtr->Get_Rxn_solution_map()) {
     // A key less than zero indicates an internal solution
@@ -47,16 +59,50 @@ void IPhreeqcPOET::parseInit() {
       continue;
     }
 
-    const essential_names curr_conc_names = this->dump_essential_names(id);
+    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<std::string> 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<std::string> solution_order(val[POET_SOL].begin() + 3,
+                                            val[POET_SOL].end());
 
     auto pair = std::make_pair(
-        curr_conc_names, this->get_essential_values(id, curr_conc_names[0]));
+        val, this->get_essential_values_init(id, solution_order));
     raw_initials[id] = pair;
   }
 }
 
-IPhreeqcPOET::essential_names IPhreeqcPOET::union_raws(const RawMap &raws) {
+std::string PhreeqcInit::subExchangeName(std::string name) {
-  IPhreeqcPOET::essential_names names;
+  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]);
@@ -71,7 +117,7 @@ IPhreeqcPOET::essential_names IPhreeqcPOET::union_raws(const RawMap &raws) {
 }
 
 std::vector<std::vector<double>>
-IPhreeqcPOET::conc_from_essentials(const RawMap &raws,
+PhreeqcInit::conc_from_essentials(const RawMap &raws,
                                   const essential_names &names) {
   std::vector<std::vector<double>> values;
 
@@ -94,10 +140,10 @@ IPhreeqcPOET::conc_from_essentials(const RawMap &raws,
   return values;
 }
 
-IPhreeqcPOET::PhreeqcMat IPhreeqcPOET::getPhreeqcMat() {
+PhreeqcInit::PhreeqcMat PhreeqcInit::buildPhreeqcMat() {
   std::vector<std::vector<double>> values;
 
-  const IPhreeqcPOET::essential_names ess_names =
+  const PhreeqcInit::essential_names ess_names =
       this->union_raws(this->raw_initials);
 
   const std::vector<std::vector<double>> conc_values =
@@ -117,8 +163,8 @@ IPhreeqcPOET::PhreeqcMat IPhreeqcPOET::getPhreeqcMat() {
   return {.names = conc_names, .ids = ids, .values = conc_values};
 }
 
-IPhreeqcPOET::ModulesArray
+PhreeqcInit::ModulesArray
-IPhreeqcPOET::getModuleSizes(const std::vector<int> &cell_ids) {
+PhreeqcInit::getModuleSizes(const std::vector<int> &cell_ids) {
   ModulesArray module_sizes;
   RawMap raws;
 
@@ -126,7 +172,7 @@ IPhreeqcPOET::getModuleSizes(const std::vector<int> &cell_ids) {
     raws[id] = this->raw_initials[id];
   }
 
-  const IPhreeqcPOET::essential_names ess_names = this->union_raws(raws);
+  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(); });

poet/src/GetSet.cpp

@@ -1,119 +0,0 @@
-#include <IPhreeqcPOET.hpp>
-
-IPhreeqcPOET::essential_names
-IPhreeqcPOET::dump_essential_names(std::size_t cell_number) {
-
-  essential_names eNames;
-
-  // Solutions
-  if (this->Get_solution(cell_number) != NULL) {
-    std::vector<std::string> &eSolNames = eNames[POET_SOL];
-    this->Get_solution(cell_number)->dump_essential_names(eSolNames);
-  }
-
-  // Exchange
-  if (this->Get_exchange(cell_number) != NULL) {
-    std::vector<std::string> &eExchNames = eNames[POET_EXCH];
-    this->Get_exchange(cell_number)->dump_essential_names(eExchNames);
-  }
-
-  // Kinetics
-  if (this->Get_kinetic(cell_number) != NULL) {
-    std::vector<std::string> &eKinNames = eNames[POET_KIN];
-    this->Get_kinetic(cell_number)->dump_essential_names(eKinNames);
-  }
-
-  // PPassemblage
-  if (this->Get_equilibrium(cell_number) != NULL) {
-    std::vector<std::string> &eEquNames = eNames[POET_EQUIL];
-    this->Get_equilibrium(cell_number)->dump_essential_names(eEquNames);
-  }
-
-  // Surface
-  if (this->Get_surface(cell_number) != NULL) {
-    std::vector<std::string> &eSurfNames = eNames[POET_SURF];
-    this->Get_surface(cell_number)->dump_essential_names(eSurfNames);
-  }
-
-  return eNames;
-}
-
-std::vector<double>
-IPhreeqcPOET::get_essential_values(std::size_t cell_number,
-                                   const std::vector<std::string> &order) {
-  std::vector<double> essentials;
-
-  // Solutions
-  if (this->Get_solution(cell_number) != NULL) {
-    std::vector<double> sol_values;
-
-    this->Get_solution(cell_number)->get_essential_values(sol_values, order);
-    essentials.insert(essentials.end(), sol_values.begin(), sol_values.end());
-  }
-
-  // Exchange
-  if (this->Get_exchange(cell_number) != NULL) {
-    std::vector<double> exch_values;
-
-    this->Get_exchange(cell_number)->get_essential_values(exch_values);
-    essentials.insert(essentials.end(), exch_values.begin(), exch_values.end());
-  }
-
-  // Kinetics
-  if (this->Get_kinetic(cell_number) != NULL) {
-    std::vector<double> kin_values;
-
-    this->Get_kinetic(cell_number)->get_essential_values(kin_values);
-    essentials.insert(essentials.end(), kin_values.begin(), kin_values.end());
-  }
-
-  // PPassemblage
-  if (this->Get_equilibrium(cell_number) != NULL) {
-    std::vector<double> equ_values;
-
-    this->Get_equilibrium(cell_number)->get_essential_values(equ_values);
-    essentials.insert(essentials.end(), equ_values.begin(), equ_values.end());
-  }
-
-  // Surface
-  if (this->Get_surface(cell_number) != NULL) {
-    std::vector<double> surf_values;
-
-    this->Get_surface(cell_number)->get_essential_values(surf_values);
-    essentials.insert(essentials.end(), surf_values.begin(), surf_values.end());
-  }
-
-  return essentials;
-}
-
-void IPhreeqcPOET::set_essential_values(std::size_t cell_number,
-                                        const std::vector<std::string> &order,
-                                        std::vector<double> &values) {
-
-  auto dump_it = values.begin();
-
-  // Solutions
-  if (this->Get_solution(cell_number) != NULL) {
-    this->Get_solution(cell_number)->set_essential_values(dump_it, order);
-  }
-
-  // Exchange
-  if (this->Get_exchange(cell_number) != NULL) {
-    this->Get_exchange(cell_number)->set_essential_values(dump_it);
-  }
-
-  // Kinetics
-  if (this->Get_kinetic(cell_number) != NULL) {
-    this->Get_kinetic(cell_number)->set_essential_values(dump_it);
-  }
-
-  // PPassemblage
-  if (this->Get_equilibrium(cell_number) != NULL) {
-    this->Get_equilibrium(cell_number)->set_essential_values(dump_it);
-  }
-
-  // Surface
-  if (this->Get_surface(cell_number) != NULL) {
-    this->Get_surface(cell_number)->set_essential_values(dump_it);
-  }
-}

poet/src/InitGetSet.cpp

@@ -0,0 +1,251 @@
+#include "EquilibriumWrapper.hpp"
+#include "ExchangeWrapper.hpp"
+#include "KineticWrapper.hpp"
+#include "PhreeqcInit.hpp"
+#include "SolutionWrapper.hpp"
+#include "SurfaceWrapper.hpp"
+#include "global_structures.h"
+
+#include <cassert>
+#include <cstddef>
+#include <set>
+#include <span>
+#include <string>
+#include <vector>
+
+std::vector<std::string> PhreeqcInit::dump_solution_names(int cell_number) {
+  constexpr std::size_t SKIP_H_O_CB = 3;
+
+  std::vector<std::string> 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<std::string> union_sol_names,
+    essential_names &names) {
+  // Exchange
+  if (this->Get_exchange(cell_number) != NULL) {
+    auto *exchange = this->Get_exchange(cell_number);
+
+    std::vector<std::string> 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<std::string> 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<std::string> 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()) {
+      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<std::string> comp_formulas;
+    std::vector<std::string> 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<std::string> PhreeqcInit::find_all_valence_states(
+    const std::vector<std::string> &solution_names, const std::size_t offset) {
+  std::vector<std::string> solution_with_valences(
+      solution_names.begin(), solution_names.begin() + offset);
+
+  // to avoid duplicates store already evaluated master species
+  std::set<std::string> 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<double>
+PhreeqcInit::get_essential_values_init(std::size_t cell_number,
+                                       const std::vector<std::string> &order) {
+  std::vector<double> essentials;
+
+  // Solutions
+  if (this->Get_solution(cell_number) != NULL) {
+    SolutionWrapper solution(this->Get_solution(cell_number), order);
+    std::vector<double> sol_values(solution.size());
+    std::span<double> 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<double> exch_values(exchange.size());
+    std::span<double> 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<double> kin_values(kinetic.size());
+    std::span<double> 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<double> equ_values(equilibrium.size());
+    std::span<double> 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<double> surf_values(surface.size());
+    std::span<double> surf_span{surf_values};
+
+    surface.get(surf_span);
+
+    essentials.insert(essentials.end(), surf_values.begin(), surf_values.end());
+  }
+
+  return essentials;
+}
+
+std::map<int, std::string> PhreeqcInit::raw_dumps() {
+  std::map<int, std::string> 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;
+}

poet/src/KineticCompWrapper.cpp

@@ -0,0 +1,40 @@
+#include "KineticWrapper.hpp"
+#include <cstddef>
+#include <string>
+#include <vector>
+
+KineticWrapper::KineticCompWrapper::KineticCompWrapper(cxxKineticsComp &comp)
+    : kin_comp(comp) {
+  this->num_elements = kin_comp.Get_d_params().size() + 1;
+}
+
+void KineticWrapper::KineticCompWrapper::get(std::span<LDBLE> &data) const {
+  data[0] = this->kin_comp.Get_m();
+
+  std::size_t next_index = 1;
+  for (const auto &param : this->kin_comp.Get_d_params()) {
+    data[next_index++] = param;
+  }
+}
+void KineticWrapper::KineticCompWrapper::set(const std::span<LDBLE> &data) {
+  this->kin_comp.Set_m(data[0]);
+
+  std::size_t next_index = 1;
+  for (auto &param : this->kin_comp.Get_d_params()) {
+    param = data[next_index++];
+  }
+}
+
+std::vector<std::string>
+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);
+
+  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/KineticWrapper.cpp

@@ -0,0 +1,78 @@
+#include "KineticWrapper.hpp"
+#include <vector>
+
+KineticWrapper::KineticWrapper(cxxKinetics *kinetics,
+                               const std::vector<std::string> &kin_comps)
+    : kinetics(kinetics) {
+
+  for (const auto &kin_name : kin_comps) {
+    auto it = std::find_if(kinetics->Get_kinetics_comps().begin(),
+                           kinetics->Get_kinetics_comps().end(),
+                           [&](const cxxKineticsComp &comp) {
+                             return comp.Get_rate_name() == kin_name;
+                           });
+
+    if (it == kinetics->Get_kinetics_comps().end()) {
+      throw std::runtime_error(
+          "Kinetic component not found in Phreeqc variables");
+    }
+
+    kinetic_comps.push_back(std::make_unique<KineticCompWrapper>(*it));
+
+    num_elements += kinetic_comps.back()->size();
+  }
+}
+
+void KineticWrapper::get(std::span<LDBLE> &data) const {
+  std::size_t offset = 0;
+
+  for (const auto &comp : kinetic_comps) {
+    std::span<LDBLE> comp_span = data.subspan(offset, comp->size());
+
+    comp->get(comp_span);
+
+    offset += comp->size();
+  }
+}
+
+void KineticWrapper::set(const std::span<LDBLE> &data) {
+  std::size_t offset = 0;
+
+  for (const auto &comp : kinetic_comps) {
+    std::span<LDBLE> comp_span = data.subspan(offset, comp->size());
+
+    comp->set(comp_span);
+
+    offset += comp->size();
+  }
+}
+
+std::vector<std::string>
+KineticWrapper::names(cxxKinetics *kinetics,
+                      std::vector<std::string> &kin_comps) {
+  std::vector<std::string> names;
+
+  for (const auto &comp : kinetics->Get_kinetics_comps()) {
+    std::vector<std::string> comp_names = KineticCompWrapper::names(comp);
+    names.insert(names.end(), comp_names.begin(), comp_names.end());
+
+    kin_comps.push_back(comp.Get_rate_name());
+  }
+
+  for (const auto &kin_name : kin_comps) {
+    auto it = std::find_if(kinetics->Get_kinetics_comps().begin(),
+                           kinetics->Get_kinetics_comps().end(),
+                           [&](const cxxKineticsComp &comp) {
+                             return comp.Get_rate_name() == kin_name;
+                           });
+
+    if (it == kinetics->Get_kinetics_comps().end()) {
+      throw std::runtime_error(
+          "Kinetic component not found in Phreeqc variables");
+    }
+
+    names.push_back(it->Get_rate_name());
+  }
+
+  return names;
+}

poet/src/SolutionWrapper.cpp

@@ -0,0 +1,63 @@
+#include "SolutionWrapper.hpp"
+#include "NameDouble.h"
+#include <vector>
+
+SolutionWrapper::SolutionWrapper(
+    cxxSolution *soln, const std::vector<std::string> &solution_order_)
+    : solution(soln), solution_order(solution_order_) {
+  this->num_elements = solution_order.size() + NUM_ESSENTIALS;
+
+  auto &totals = solution->Get_totals();
+}
+
+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();
+
+  std::size_t i = NUM_ESSENTIALS;
+  for (const auto &tot_name : solution_order) {
+    auto it = solution->Get_totals().find(tot_name);
+    if (it == solution->Get_totals().end()) {
+      data[i++] = 0.0;
+      continue;
+    }
+    data[i++] = it->second;
+  }
+}
+
+void SolutionWrapper::set(const std::span<LDBLE> &data) {
+  std::size_t i = NUM_ESSENTIALS;
+  cxxNameDouble new_totals;
+  for (const auto &tot_name : solution_order) {
+    const double value = data[i++];
+
+    if (value < 1E-25) {
+      continue;
+    }
+    new_totals[tot_name] = value;
+  }
+
+  this->solution->Update(data[0], data[1], data[2], new_totals);
+}
+
+std::vector<std::string>
+SolutionWrapper::names(cxxSolution *solution,
+                       std::vector<std::string> &solution_order) {
+  std::vector<std::string> names;
+
+  names.push_back("H");
+  names.push_back("O");
+  names.push_back("Charge");
+
+  for (const auto &[tot_name, value] : solution->Get_totals()) {
+    if (tot_name == "H(0)" || tot_name == "O(0)") {
+      continue;
+    }
+
+    solution_order.push_back(tot_name);
+  }
+
+  names.insert(names.end(), solution_order.begin(), solution_order.end());
+  return names;
+}

poet/src/SurfaceChargeWrapper.cpp

@@ -0,0 +1,72 @@
+#include "SurfaceWrapper.hpp"
+
+SurfaceWrapper::SurfaceChargeWrapper::SurfaceChargeWrapper(
+    cxxSurfaceCharge &charge, const std::set<std::string> &solution_primaries)
+    : surface_charge(charge), primaries(solution_primaries) {
+  num_elements = NUM_NOT_TOTALS + solution_primaries.size();
+}
+
+void SurfaceWrapper::SurfaceChargeWrapper::get(
+    std::span<LDBLE> &surface) const {
+  surface[0] = this->surface_charge.Get_specific_area();
+  surface[1] = this->surface_charge.Get_grams();
+  surface[2] = this->surface_charge.Get_charge_balance();
+  surface[3] = this->surface_charge.Get_mass_water();
+  surface[4] = this->surface_charge.Get_la_psi();
+
+  const auto dl_map = this->surface_charge.Get_diffuse_layer_totals();
+
+  for (auto it = this->primaries.begin(); it != this->primaries.end(); it++) {
+    const std::size_t index =
+        NUM_NOT_TOTALS + std::distance(this->primaries.begin(), it);
+
+    const auto dl_it = dl_map.find(*it);
+
+    if (dl_it == dl_map.end()) {
+      surface[index] = 0;
+      continue;
+    }
+    surface[index] = dl_it->second;
+  }
+}
+
+void SurfaceWrapper::SurfaceChargeWrapper::set(
+    const std::span<LDBLE> &surface) {
+  this->surface_charge.Set_specific_area(surface[0]);
+  this->surface_charge.Set_grams(surface[1]);
+  this->surface_charge.Set_charge_balance(surface[2]);
+  this->surface_charge.Set_mass_water(surface[3]);
+  this->surface_charge.Set_la_psi(surface[4]);
+
+  auto &dl_map = this->surface_charge.Get_diffuse_layer_totals();
+  dl_map.clear();
+
+  for (auto it = this->primaries.begin(); it != this->primaries.end(); it++) {
+    const std::size_t index =
+        NUM_NOT_TOTALS + std::distance(this->primaries.begin(), it);
+
+    if (surface[index] == 0) {
+      continue;
+    }
+
+    dl_map[*it] = surface[index];
+  }
+}
+
+std::vector<std::string> SurfaceWrapper::SurfaceChargeWrapper::names(
+    cxxSurfaceCharge *charge, const std::set<std::string> &primaries) {
+  std::vector<std::string> names;
+
+  const std::string &charge_name = charge->Get_name();
+  names.push_back(charge_name + "_area");
+  names.push_back(charge_name + "_grams");
+  names.push_back(charge_name + "_cb");
+  names.push_back(charge_name + "_mw");
+  names.push_back(charge_name + "_la");
+
+  for (const auto &name : primaries) {
+    names.push_back(charge_name + "_tot_" + name);
+  }
+
+  return names;
+}

poet/src/SurfaceCompWrapper.cpp

@@ -0,0 +1,52 @@
+#include "SurfaceWrapper.hpp"
+
+SurfaceWrapper::SurfaceCompWrapper::SurfaceCompWrapper(cxxSurfaceComp &comp)
+    : surface_comp(comp) {
+  num_elements = NUM_NOT_TOTALS + comp.Get_totals().size();
+
+  for (const auto &[name, value] : comp.Get_totals()) {
+    total_names.push_back(name);
+  }
+}
+
+void SurfaceWrapper::SurfaceCompWrapper::get(std::span<LDBLE> &surface) const {
+  surface[0] = this->surface_comp.Get_moles();
+  surface[1] = this->surface_comp.Get_la();
+  surface[2] = this->surface_comp.Get_charge_balance();
+
+  const auto &totals = this->surface_comp.Get_totals();
+  for (std::size_t i = 0; i < this->total_names.size(); i++) {
+    surface[NUM_NOT_TOTALS + i] = totals.at(this->total_names[i]);
+  }
+}
+
+void SurfaceWrapper::SurfaceCompWrapper::set(const std::span<LDBLE> &surface) {
+  this->surface_comp.Set_moles(surface[0]);
+  this->surface_comp.Set_la(surface[1]);
+  this->surface_comp.Set_charge_balance(surface[2]);
+
+  auto &totals = this->surface_comp.Get_totals();
+  totals.clear();
+
+  for (std::size_t i = 0; i < this->total_names.size(); i++) {
+    const std::size_t index = NUM_NOT_TOTALS + i;
+
+    totals[this->total_names[i]] = surface[index];
+  }
+}
+
+std::vector<std::string>
+SurfaceWrapper::SurfaceCompWrapper::names(cxxSurfaceComp &comp) {
+  std::vector<std::string> names;
+
+  const std::string &phase_name = comp.Get_formula();
+  names.push_back(phase_name + "_moles");
+  names.push_back(phase_name + "_la");
+  names.push_back(phase_name + "_cb");
+
+  for (const auto &tot : comp.Get_totals()) {
+    names.push_back(phase_name + "_" + tot.first);
+  }
+
+  return names;
+}

poet/src/SurfaceWrapper.cpp

@@ -0,0 +1,98 @@
+#include "SurfaceWrapper.hpp"
+#include "SurfaceComp.h"
+
+SurfaceWrapper::SurfaceWrapper(cxxSurface *surf,
+                               const std::set<std::string> &solution_primaries,
+                               const std::vector<std::string> &comp_formulas,
+                               const std::vector<std::string> &charge_names)
+    : surface(surf) {
+
+  for (const auto &comp_name : comp_formulas) {
+    auto it = std::find_if(surface->Get_surface_comps().begin(),
+                           surface->Get_surface_comps().end(),
+                           [&](const cxxSurfaceComp &comp) {
+                             return comp.Get_formula() == comp_name;
+                           });
+
+    if (it == surface->Get_surface_comps().end()) {
+      throw std::runtime_error(
+          "Surface component not found in Phreeqc variables");
+    }
+
+    surface_comps.push_back(std::make_unique<SurfaceCompWrapper>(*it));
+
+    num_elements += surface_comps.back()->size();
+  }
+
+  for (const auto &charge_name : charge_names) {
+    auto it = std::find_if(surface->Get_surface_charges().begin(),
+                           surface->Get_surface_charges().end(),
+                           [&](const cxxSurfaceCharge &charge) {
+                             return charge.Get_name() == charge_name;
+                           });
+
+    if (it == surface->Get_surface_charges().end()) {
+      throw std::runtime_error("Surface charge not found in Phreeqc variables");
+    }
+
+    surface_charges.push_back(
+        std::make_unique<SurfaceChargeWrapper>(*it, solution_primaries));
+
+    num_elements += surface_charges.back()->size();
+  }
+}
+
+void SurfaceWrapper::get(std::span<LDBLE> &surface) const {
+  std::size_t offset = 0;
+
+  for (const auto &comp : surface_comps) {
+    std::span<LDBLE> comp_span = surface.subspan(offset, comp->size());
+    comp->get(comp_span);
+    offset += comp->size();
+  }
+
+  for (const auto &charge : surface_charges) {
+    std::span<LDBLE> charge_span = surface.subspan(offset, charge->size());
+    charge->get(charge_span);
+    offset += charge->size();
+  }
+}
+
+void SurfaceWrapper::set(const std::span<LDBLE> &surface) {
+  std::size_t offset = 0;
+
+  for (const auto &comp : surface_comps) {
+    std::span<LDBLE> comp_span = surface.subspan(offset, comp->size());
+    comp->set(comp_span);
+    offset += comp->size();
+  }
+
+  for (const auto &charge : surface_charges) {
+    std::span<LDBLE> charge_span = surface.subspan(offset, charge->size());
+    charge->set(charge_span);
+    offset += charge->size();
+  }
+}
+
+std::vector<std::string>
+SurfaceWrapper::names(cxxSurface *surface,
+                      const std::set<std::string> &solution_primaries,
+                      std::vector<std::string> &comp_formulas,
+                      std::vector<std::string> &charge_names) {
+  std::vector<std::string> names;
+
+  for (auto &comp : surface->Get_surface_comps()) {
+    comp_formulas.push_back(comp.Get_formula());
+    const auto charge_names = SurfaceCompWrapper::names(comp);
+    names.insert(names.end(), charge_names.begin(), charge_names.end());
+  }
+
+  for (auto &charge : surface->Get_surface_charges()) {
+    charge_names.push_back(charge.Get_name());
+    const auto comp_names =
+        SurfaceChargeWrapper::names(&charge, solution_primaries);
+    names.insert(names.end(), comp_names.begin(), comp_names.end());
+  }
+
+  return names;
+}

src/phreeqcpp/Exchange.cxx

@@ -2,17 +2,18 @@
 //
 //////////////////////////////////////////////////////////////////////
 #ifdef _DEBUG
-#pragma warning(disable : 4786)	// disable truncation warning (Only used by debugger)
+#pragma warning(                                                               \
+    disable : 4786) // disable truncation warning (Only used by debugger)
 #endif
-#include <iostream>				// std::cout std::cerr
-#include <cassert>				// assert
 #include <algorithm> // std::sort
+#include <cassert>   // assert
+#include <iostream>  // std::cout std::cerr
 
 #include "Utils.h" // define first
 
+#include "Exchange.h"
 #include "Phreeqc.h"
 #include "cxxMix.h"
-#include "Exchange.h"
 #include "phqalloc.h"
 
 #if defined(PHREEQCI_GUI)
@@ -31,72 +32,54 @@ cxxExchange::cxxExchange(PHRQ_io *io)
     //
     // default constructor for cxxExchange
     //
-:	cxxNumKeyword(io)
+    : cxxNumKeyword(io) {
-{
   new_def = false;
   solution_equilibria = false;
   n_solution = -999;
   pitzer_exchange_gammas = true;
 }
-cxxExchange::cxxExchange(const std::map < int, cxxExchange > &entities,
+cxxExchange::cxxExchange(const std::map<int, cxxExchange> &entities,
-						 cxxMix & mix, int l_n_user, PHRQ_io *io):
+                         cxxMix &mix, int l_n_user, PHRQ_io *io)
-cxxNumKeyword(io)
+    : cxxNumKeyword(io) {
-{
   this->n_user = this->n_user_end = l_n_user;
   this->pitzer_exchange_gammas = true;
   this->new_def = false;
   this->n_solution = -999;
   this->solution_equilibria = false;
-//
+  //
-//   Mix exchangers
+  //   Mix exchangers
-//
+  //
-	const std::map < int, LDBLE >&mixcomps = mix.Get_mixComps();
+  const std::map<int, LDBLE> &mixcomps = mix.Get_mixComps();
-	std::map < int, LDBLE >::const_iterator it;
+  std::map<int, LDBLE>::const_iterator it;
-	for (it = mixcomps.begin(); it != mixcomps.end(); it++)
+  for (it = mixcomps.begin(); it != mixcomps.end(); it++) {
-	{
+    if (entities.find(it->first) != entities.end()) {
-		if (entities.find(it->first) != entities.end())
+      const cxxExchange *entity_ptr = &(entities.find(it->first)->second);
-		{
-			const cxxExchange *entity_ptr =
-				&(entities.find(it->first)->second);
       this->add(*entity_ptr, it->second);
       this->pitzer_exchange_gammas = entity_ptr->pitzer_exchange_gammas;
     }
   }
 }
 
-
+cxxExchange::~cxxExchange() {}
-cxxExchange::~cxxExchange()
+bool cxxExchange::Get_related_phases() const {
-{
+  for (size_t i = 0; i < this->exchange_comps.size(); i++) {
-}
+    if (this->exchange_comps[i].Get_phase_name().size() > 0) {
-bool
-cxxExchange::Get_related_phases() const
-{
-	for (size_t i = 0; i < this->exchange_comps.size(); i++)
-	{
-		if (this->exchange_comps[i].Get_phase_name().size() > 0)
-		{
       return (true);
     }
   }
   return (false);
 }
 
-bool
+bool cxxExchange::Get_related_rate() const {
-cxxExchange::Get_related_rate() const
+  for (size_t i = 0; i < this->exchange_comps.size(); i++) {
-{
+    if (this->exchange_comps[i].Get_rate_name().size() > 0) {
-	for (size_t i = 0; i < this->exchange_comps.size(); i++)
-	{
-		if (this->exchange_comps[i].Get_rate_name().size() > 0)
-		{
       return (true);
     }
   }
   return (false);
 }
 
-void
+void cxxExchange::dump_xml(std::ostream &s_oss, unsigned int indent) const {
-cxxExchange::dump_xml(std::ostream & s_oss, unsigned int indent) const
-{
   unsigned int i;
   s_oss.precision(DBL_DIG - 1);
   std::string indent0(""), indent1(""), indent2("");
@@ -109,25 +92,25 @@ cxxExchange::dump_xml(std::ostream & s_oss, unsigned int indent) const
 
   // Exchange element and attributes
   s_oss << indent0;
-	s_oss << "<exchange " << "\n";
+  s_oss << "<exchange "
+        << "\n";
 
   s_oss << indent1;
-	s_oss << "pitzer_exchange_gammas=\"" << this->
+  s_oss << "pitzer_exchange_gammas=\"" << this->pitzer_exchange_gammas << "\""
-		pitzer_exchange_gammas << "\"" << "\n";
+        << "\n";
 
   // components
   s_oss << indent1;
-	s_oss << "<component " << "\n";
+  s_oss << "<component "
-	for (size_t j = 0; j < this->exchange_comps.size(); j++)
+        << "\n";
-	{
+  for (size_t j = 0; j < this->exchange_comps.size(); j++) {
     this->exchange_comps[j].dump_xml(s_oss, indent + 2);
   }
 
   return;
 }
-void
+void cxxExchange::dump_raw(std::ostream &s_oss, unsigned int indent,
-cxxExchange::dump_raw(std::ostream & s_oss, unsigned int indent, int *n_out) const
+                           int *n_out) const {
-{
   unsigned int i;
   s_oss.precision(DBL_DIG - 1);
   std::string indent0(""), indent1(""), indent2("");
@@ -141,16 +124,18 @@ cxxExchange::dump_raw(std::ostream & s_oss, unsigned int indent, int *n_out) con
   // Exchange element and attributes
   s_oss << indent0;
   int n_user_local = (n_out != NULL) ? *n_out : this->n_user;
-	s_oss << "EXCHANGE_RAW                 " << n_user_local << " " << this->description << "\n";
+  s_oss << "EXCHANGE_RAW                 " << n_user_local << " "
+        << this->description << "\n";
 
   s_oss << indent1 << "# EXCHANGE_MODIFY candidate identifiers #\n";
   s_oss << indent1;
-	s_oss << "-exchange_gammas           " << (this->pitzer_exchange_gammas ? 1 : 0) << "\n";
+  s_oss << "-exchange_gammas           "
+        << (this->pitzer_exchange_gammas ? 1 : 0) << "\n";
   // exchComps
-	for (size_t i = 0; i < this->exchange_comps.size(); i++)
+  for (size_t i = 0; i < this->exchange_comps.size(); i++) {
-	{
     s_oss << indent1;
-		s_oss << "-component                 " << this->exchange_comps[i].Get_formula() << "\n";
+    s_oss << "-component                 "
+          << this->exchange_comps[i].Get_formula() << "\n";
     this->exchange_comps[i].dump_raw(s_oss, indent + 2);
   }
 
@@ -164,68 +149,14 @@ cxxExchange::dump_raw(std::ostream & s_oss, unsigned int indent, int *n_out) con
 
   s_oss << indent1 << "# Exchange workspace variables #\n";
   s_oss << indent1;
-	s_oss << "-totals" << "\n";
+  s_oss << "-totals"
+        << "\n";
   this->totals.dump_raw(s_oss, indent + 1);
 
   return;
 }
 
-void cxxExchange::dump_essential_names(
+void cxxExchange::read_raw(CParser &parser, bool check) {
-    std::vector<std::string> &e_names) const {
-  e_names.clear();
-  e_names.reserve(this->exchange_comps.size() *
-                  this->exchange_comps.at(0).Get_totals().size());
-
-  for (const auto &comp : this->exchange_comps) {
-    const std::string formular = comp.Get_formula();
-    e_names.push_back(formular);
-    for (const auto &total : comp.Get_totals()) {
-      if (total.first == formular) {
-        continue;
-      }
-      e_names.push_back(total.first + formular);
-    }
-  }
-}
-
-void cxxExchange::get_essential_values(std::vector<LDBLE> &e_values) const {
-  e_values.clear();
-  e_values.reserve(this->exchange_comps.size() *
-                   this->exchange_comps.at(0).Get_totals().size());
-
-  for (const auto &comp : this->exchange_comps) {
-    const std::string formular = comp.Get_formula();
-    e_values.push_back(comp.Get_totals().find(formular)->second);
-    for (const auto &total : comp.Get_totals()) {
-      if (total.first == formular) {
-        continue;
-      }
-      e_values.push_back(total.second);
-    }
-  }
-
-  return;
-}
-
-void cxxExchange::set_essential_values(std::vector<LDBLE>::iterator &it) {
-  for (auto &comp : this->exchange_comps) {
-    const std::string formular = comp.Get_formula();
-    auto formular_total = comp.Get_totals();
-    formular_total[formular] = *(it++);
-    for (auto &total : comp.Get_totals()) {
-      if (total.first == formular) {
-        continue;
-      }
-      total.second = *(it++);
-    }
-  }
-
-  return;
-}
-
-void
-cxxExchange::read_raw(CParser & parser, bool check)
-{
   std::istream::pos_type ptr;
   std::istream::pos_type next_char;
   std::string token;
@@ -237,20 +168,15 @@ cxxExchange::read_raw(CParser & parser, bool check)
 
   bool pitzer_exchange_gammas_defined(false);
 
-	for (;;)
+  for (;;) {
-	{
     int opt;
-		if (useLastLine == false)
+    if (useLastLine == false) {
-		{
       opt = parser.get_option(vopts, next_char);
-		}
+    } else {
-		else
-		{
       opt = parser.getOptionFromLastLine(vopts, next_char, true);
     }
     useLastLine = false;
-		switch (opt)
+    switch (opt) {
-		{
     case CParser::OPT_EOF:
       break;
     case CParser::OPT_KEYWORD:
@@ -265,13 +191,10 @@ cxxExchange::read_raw(CParser & parser, bool check)
 
     case 0: // pitzer_exchange_gammas
     case 2: // exchange_gammas
-			if (!(parser.get_iss() >> this->pitzer_exchange_gammas))
+      if (!(parser.get_iss() >> this->pitzer_exchange_gammas)) {
-			{
         this->pitzer_exchange_gammas = false;
         parser.incr_input_error();
-				parser.
+        parser.error_msg("Expected boolean value for pitzer_exchange_gammas.",
-					error_msg
-					("Expected boolean value for pitzer_exchange_gammas.",
                          PHRQ_io::OT_CONTINUE);
       }
       pitzer_exchange_gammas_defined = true;
@@ -279,81 +202,61 @@ cxxExchange::read_raw(CParser & parser, bool check)
     case 1: // component
     {
       std::string str;
-				if (!(parser.get_iss() >> str))
+      if (!(parser.get_iss() >> str)) {
-				{
         parser.incr_input_error();
         parser.error_msg("Expected string value for component name.",
                          PHRQ_io::OT_CONTINUE);
-				}
+      } else {
-				else
-				{
         cxxExchComp temp_comp(this->io);
         temp_comp.Set_formula(str.c_str());
         cxxExchComp *comp_ptr = this->Find_comp(str);
-					if (comp_ptr)
+        if (comp_ptr) {
-					{
           temp_comp = *comp_ptr;
         }
         temp_comp.read_raw(parser, check);
-					if (comp_ptr)
+        if (comp_ptr) {
-					{
+          for (size_t j = 0; j < this->exchange_comps.size(); j++) {
-						for (size_t j = 0; j < this->exchange_comps.size(); j++)
+            if (Utilities::strcmp_nocase(
-						{
+                    this->exchange_comps[j].Get_formula().c_str(),
-							if (Utilities::strcmp_nocase(this->exchange_comps[j].Get_formula().c_str(), str.c_str()) == 0)
+                    str.c_str()) == 0) {
-							{
               this->exchange_comps[j] = temp_comp;
             }
           }
-					}
+        } else {
-					else
-					{
           this->exchange_comps.push_back(temp_comp);
         }
         useLastLine = true;
       }
-			}
+    } break;
-			break;
     case 3: // new_def
-			if (!(parser.get_iss() >> this->new_def))
+      if (!(parser.get_iss() >> this->new_def)) {
-			{
         this->new_def = false;
         parser.incr_input_error();
-				parser.
+        parser.error_msg("Expected boolean value for new_def.",
-					error_msg
-					("Expected boolean value for new_def.",
                          PHRQ_io::OT_CONTINUE);
       }
       break;
     case 4: // solution_equilibria
-			if (!(parser.get_iss() >> this->solution_equilibria))
+      if (!(parser.get_iss() >> this->solution_equilibria)) {
-			{
         this->solution_equilibria = false;
         parser.incr_input_error();
-				parser.
+        parser.error_msg("Expected boolean value for solution_equilibria.",
-					error_msg
-					("Expected boolean value for solution_equilibria.",
                          PHRQ_io::OT_CONTINUE);
       }
       break;
     case 5: // n_solution
-			if (!(parser.get_iss() >> this->n_solution))
+      if (!(parser.get_iss() >> this->n_solution)) {
-			{
         this->n_solution = -999;
         parser.incr_input_error();
-				parser.
+        parser.error_msg("Expected integer value for n_solution.",
-					error_msg
-					("Expected integer value for n_solution.",
                          PHRQ_io::OT_CONTINUE);
       }
       break;
     case 6: // totals
-			if (this->totals.read_raw(parser, next_char) !=
+      if (this->totals.read_raw(parser, next_char) != CParser::PARSER_OK) {
-				CParser::PARSER_OK)
-			{
         parser.incr_input_error();
-				parser.
+        parser.error_msg(
-					error_msg
+            "Expected element name and molality for Exchange totals.",
-					("Expected element name and molality for Exchange totals.",
             PHRQ_io::OT_CONTINUE);
       }
       break;
@@ -362,39 +265,33 @@ cxxExchange::read_raw(CParser & parser, bool check)
       break;
   }
   // members that must be defined
-	if (check)
+  if (check) {
-	{
+    if (pitzer_exchange_gammas_defined == false) {
-		if (pitzer_exchange_gammas_defined == false)
-		{
       parser.incr_input_error();
-			parser.
+      parser.error_msg(
-				error_msg
+          "Pitzer_exchange_gammsa not defined for EXCHANGE_RAW input.",
-				("Pitzer_exchange_gammsa not defined for EXCHANGE_RAW input.",
           PHRQ_io::OT_CONTINUE);
     }
   }
   this->Sort_comps();
 }
-void
+void cxxExchange::add(const cxxExchange &addee, LDBLE extensive)
-cxxExchange::add(const cxxExchange & addee, LDBLE extensive)
+//
-		//
+// Add existing exchange to "this" exchange
-		// Add existing exchange to "this" exchange
+//
-		//
 {
   // exchComps
   if (extensive == 0.0)
     return;
-	for (size_t i = 0; i < addee.exchange_comps.size(); i++)
+  for (size_t i = 0; i < addee.exchange_comps.size(); i++) {
-	{
     size_t j;
     for (j = 0; j < this->Get_exchange_comps().size(); j++)
-		if (addee.exchange_comps[i].Get_formula() == this->exchange_comps[j].Get_formula())
+      if (addee.exchange_comps[i].Get_formula() ==
-		{
+          this->exchange_comps[j].Get_formula()) {
         this->exchange_comps[j].add(addee.exchange_comps[i], extensive);
         break;
       }
-		if (j == this->exchange_comps.size())
+    if (j == this->exchange_comps.size()) {
-		{
       cxxExchComp exc = addee.exchange_comps[i];
       exc.multiply(extensive);
       this->exchange_comps.push_back(exc);
@@ -402,76 +299,56 @@ cxxExchange::add(const cxxExchange & addee, LDBLE extensive)
   }
   this->pitzer_exchange_gammas = addee.pitzer_exchange_gammas;
 }
-void
+void cxxExchange::totalize() {
-cxxExchange::totalize()
-{
   this->totals.clear();
   // component structures
-	for (size_t i = 0; i < this->exchange_comps.size(); i++)
+  for (size_t i = 0; i < this->exchange_comps.size(); i++) {
-	{
     this->totals.add_extensive(this->exchange_comps[i].Get_totals(), 1.0);
     this->totals.add("Charge", this->exchange_comps[i].Get_charge_balance());
   }
   return;
 }
-bool 
+bool cxxExchange::Get_pitzer_exchange_gammas() const {
-cxxExchange::Get_pitzer_exchange_gammas() const
-{
   return this->pitzer_exchange_gammas;
 }
-void 
+void cxxExchange::Set_pitzer_exchange_gammas(bool b) {
-cxxExchange::Set_pitzer_exchange_gammas(bool b)
-{
   this->pitzer_exchange_gammas = b;
 }
-const cxxNameDouble & 
+const cxxNameDouble &cxxExchange::Get_totals() const { return totals; }
-cxxExchange::Get_totals() const
+cxxExchComp *cxxExchange::Find_comp(std::string s) {
-{
+  for (size_t i = 0; i < this->exchange_comps.size(); i++) {
-	return totals;
-}
-cxxExchComp *cxxExchange::Find_comp(std::string s)
-{
-	for (size_t i = 0; i < this->exchange_comps.size(); i++)
-	{
     cxxNameDouble nd(this->exchange_comps[i].Get_totals());
     cxxNameDouble::iterator nd_it;
-		for (nd_it = nd.begin(); nd_it != nd.end(); nd_it++)
+    for (nd_it = nd.begin(); nd_it != nd.end(); nd_it++) {
-		{
+      if (nd_it->first == s) {
-			if(nd_it->first == s)
-			{
         return (&(this->exchange_comps[i]));
       }
     }
   }
   return NULL;
 }
-void cxxExchange::
+void cxxExchange::Sort_comps(void) {
-Sort_comps(void)
-{
   // sort comps
   {
     std::map<std::string, cxxExchComp> comp_map;
-		for (size_t i = 0; i < this->exchange_comps.size(); i++)
+    for (size_t i = 0; i < this->exchange_comps.size(); i++) {
-		{
       comp_map[this->exchange_comps[i].Get_formula()] = this->exchange_comps[i];
     }
     this->exchange_comps.clear();
     std::map<std::string, cxxExchComp>::iterator it;
-		for (it = comp_map.begin(); it != comp_map.end(); it++)
+    for (it = comp_map.begin(); it != comp_map.end(); it++) {
-		{
       this->exchange_comps.push_back(it->second);
     }
   }
 }
 /* ---------------------------------------------------------------------- */
-void
+void cxxExchange::Serialize(Dictionary &dictionary, std::vector<int> &ints,
-cxxExchange::Serialize(Dictionary & dictionary, std::vector < int >&ints, std::vector < double >&doubles)
+                            std::vector<double> &doubles)
 /* ---------------------------------------------------------------------- */
 {
   ints.push_back(this->n_user);
-	ints.push_back((int) this->exchange_comps.size());
+  ints.push_back((int)this->exchange_comps.size());
-	for (size_t i = 0; i < this->exchange_comps.size(); i++)
+  for (size_t i = 0; i < this->exchange_comps.size(); i++) {
-	{
     exchange_comps[i].Serialize(dictionary, ints, doubles);
   }
   ints.push_back(this->pitzer_exchange_gammas ? 1 : 0);
@@ -479,12 +356,11 @@ cxxExchange::Serialize(Dictionary & dictionary, std::vector < int >&ints, std::v
   ints.push_back(this->solution_equilibria ? 1 : 0);
   ints.push_back(this->n_solution);
   this->totals.Serialize(dictionary, ints, doubles);
-
 }
 
 /* ---------------------------------------------------------------------- */
-void
+void cxxExchange::Deserialize(Dictionary &dictionary, std::vector<int> &ints,
-cxxExchange::Deserialize(Dictionary & dictionary, std::vector < int >&ints, std::vector < double >&doubles, int &ii, int &dd)
+                              std::vector<double> &doubles, int &ii, int &dd)
 /* ---------------------------------------------------------------------- */
 {
   this->n_user = ints[ii++];
@@ -493,8 +369,7 @@ cxxExchange::Deserialize(Dictionary & dictionary, std::vector < int >&ints, std:
 
   int count = ints[ii++];
   this->exchange_comps.clear();
-	for (int n = 0; n < count; n++)
+  for (int n = 0; n < count; n++) {
-	{
     cxxExchComp ec(this->io);
     ec.Deserialize(dictionary, ints, doubles, ii, dd);
     this->exchange_comps.push_back(ec);
@@ -504,17 +379,17 @@ cxxExchange::Deserialize(Dictionary & dictionary, std::vector < int >&ints, std:
   this->solution_equilibria = (ints[ii++] != 0);
   this->n_solution = ints[ii++];
   this->totals.Deserialize(dictionary, ints, doubles, ii, dd);
-
 }
 
-
+const std::vector<std::string>::value_type temp_vopts[] = {
-const std::vector< std::string >::value_type temp_vopts[] = {
+    std::vector<std::string>::value_type("pitzer_exchange_gammas"), // 0
-	std::vector< std::string >::value_type("pitzer_exchange_gammas"),   // 0
+    std::vector<std::string>::value_type("component"),              // 1
-	std::vector< std::string >::value_type("component"),		        // 1
+    std::vector<std::string>::value_type("exchange_gammas"),        // 2
-	std::vector< std::string >::value_type("exchange_gammas"),	        // 2
+    std::vector<std::string>::value_type("new_def"),                // 3
-	std::vector< std::string >::value_type("new_def"),		            // 3
+    std::vector<std::string>::value_type("solution_equilibria"),    // 4
-	std::vector< std::string >::value_type("solution_equilibria"),	    // 4
+    std::vector<std::string>::value_type("n_solution"),             // 5
-	std::vector< std::string >::value_type("n_solution"),		        // 5
+    std::vector<std::string>::value_type("totals")                  // 6
-	std::vector< std::string >::value_type("totals")		            // 6
 };
-const std::vector< std::string > cxxExchange::vopts(temp_vopts, temp_vopts + sizeof temp_vopts / sizeof temp_vopts[0]);
+const std::vector<std::string>
+    cxxExchange::vopts(temp_vopts,
+                       temp_vopts + sizeof temp_vopts / sizeof temp_vopts[0]);

src/phreeqcpp/Exchange.h

@@ -20,11 +20,8 @@ public:
 		cxxMix & mx, int n_user, PHRQ_io *io=NULL);
 	~cxxExchange();
 
-	void dump_raw(std::ostream & s_oss, unsigned int indent, int *n_out=NULL) const;
+        void dump_raw(std::ostream &s_oss, unsigned int indent,
-        void dump_essential_names(std::vector<std::string> &e_names) const;
+                      int *n_out = NULL) const;
-
-        void get_essential_values(std::vector<LDBLE> &e_values) const;
-        void set_essential_values(std::vector<LDBLE>::iterator &it);
 
         void read_raw(CParser & parser, bool check = true);
 

src/phreeqcpp/PPassemblage.cxx

@@ -131,40 +131,6 @@ cxxPPassemblage::dump_raw(std::ostream & s_oss, unsigned int indent, int *n_out)
 	this->assemblage_totals.dump_raw(s_oss, indent + 1);
 }
 
-void cxxPPassemblage::dump_essential_names(
-    std::vector<std::string> &e_names) const {
-  e_names.clear();
-  e_names.reserve(this->pp_assemblage_comps.size() * 2);
-
-  for (const auto &pp_comp : this->pp_assemblage_comps) {
-    e_names.push_back(pp_comp.first);
-    e_names.push_back(pp_comp.first + "_si");
-  }
-
-  return;
-}
-
-void cxxPPassemblage::get_essential_values(std::vector<LDBLE> &e_values) const {
-  e_values.clear();
-  e_values.reserve(this->pp_assemblage_comps.size() * 2);
-
-  for (const auto &pp_comp : this->pp_assemblage_comps) {
-    e_values.push_back(pp_comp.second.Get_moles());
-    e_values.push_back(pp_comp.second.Get_si());
-  }
-
-  return;
-}
-
-void cxxPPassemblage::set_essential_values(std::vector<LDBLE>::iterator &it) {
-  for (auto &pp_comp : this->pp_assemblage_comps) {
-    pp_comp.second.Set_moles(*(it++));
-    pp_comp.second.Set_si(*(it++));
-  }
-
-  return;
-}
-
 void
 cxxPPassemblage::read_raw(CParser & parser, bool check)
 {

src/phreeqcpp/PPassemblage.h

@@ -23,10 +23,6 @@ class cxxPPassemblage:public cxxNumKeyword
 
 	void dump_raw(std::ostream & s_oss, unsigned int indent, int *n_out=NULL) const;
 
-        void dump_essential_names(std::vector<std::string> &e_names) const;
-        void get_essential_values(std::vector<LDBLE> &e_values) const;
-        void set_essential_values(std::vector<LDBLE>::iterator &it);
-
         void read_raw(CParser &parser, bool check = true);
 
         const cxxNameDouble & Get_assemblage_totals() const

src/phreeqcpp/Phreeqc.h

@@ -424,6 +424,9 @@ public:
 	int initial_exchangers(int print);
 	int initial_gas_phases(int print);
 	int initial_solutions(int print);
+
+        int initial_solutions_poet(int sol_id);
+
         int step_save_exch(int n_user);
 	int step_save_surf(int n_user);
 	int initial_surfaces(int print);
@@ -1830,6 +1833,8 @@ protected:
 	friend class IPhreeqcMMS;
 	friend class IPhreeqcPhast;
 	friend class PhreeqcRM;
+        friend class PhreeqcInit;
+        friend class PhreeqcEngine;
 
         std::vector<int> keycount;  // used to mark keywords that have been read 
 

src/phreeqcpp/Solution.cxx

@@ -367,70 +367,6 @@ cxxSolution::dump_raw(std::ostream & s_oss, unsigned int indent, int *n_out) con
 	return;
 }
 
-void cxxSolution::dump_essential_names(
-    std::vector<std::string> &e_names) const {
-  e_names.clear();
-  e_names.reserve(this->totals.size() + 3);
-
-  // e_names.push_back("H20");
-  e_names.push_back("H");
-  e_names.push_back("O");
-  e_names.push_back("Charge");
-
-  for (const auto &total : this->totals) {
-    if ((total.first == "O(0)") || (total.first == "H(0)")) {
-      continue;
-    }
-    e_names.push_back(total.first);
-  }
-
-  return;
-}
-
-void cxxSolution::get_essential_values(
-    std::vector<LDBLE> &e_values, const std::vector<std::string> &order) const {
-  e_values.clear();
-  e_values.reserve(this->totals.size() + 4);
-
-  // e_values.push_back(mass_water);
-  e_values.push_back(total_h);
-  e_values.push_back(total_o);
-  e_values.push_back(cb);
-
-  for (int i = 3; i < order.size(); i++) {
-    const auto tot = totals.find(order[i]);
-    if (tot == totals.end()) {
-      e_values.push_back(0);
-      continue;
-    }
-    e_values.push_back((tot->second));
-  }
-
-  return;
-}
-
-void cxxSolution::set_essential_values(std::vector<LDBLE>::iterator &it,
-                                       const std::vector<std::string> &order) {
-  // skip first field of moles H2O
-  //  it++;
-
-  double t_h = *(it++);
-  double t_o = *(it++);
-  double c_b = *(it++);
-
-  cxxNameDouble nd;
-  for (int i = 3; i < order.size(); i++) {
-    double curr_val = *(it++);
-
-    if (curr_val < 0) {
-      curr_val = 0;
-    }
-    nd.add(order[i].c_str(), curr_val);
-  }
-
-  this->Update(t_h, t_o, c_b, nd);
-}
-
 #ifdef USE_REVISED_READ_RAW
 void
 cxxSolution::read_raw(CParser & parser, bool check)
@@ -1259,16 +1195,19 @@ cxxSolution::Update(LDBLE h_tot, LDBLE o_tot, LDBLE charge, const cxxNameDouble
 	this->mass_water = o_tot / 55.5;
 
         // Don`t bother to update activities?
-	this->Update(const_nd);
+        // this->Update(const_nd);
-	//this->totals = const_nd;
+
-	cxxNameDouble::iterator it;
+        // UP: instead of updating, we just set the new totals ...
-	for (it = this->totals.begin(); it != this->totals.end(); it++)
+        this->totals = const_nd;
-	{
+
-		if (it->second < 1e-25)
+        // UP: ... and doesn't need another check for zero values
-		{
+
-			it->second = 0.0;
+        // cxxNameDouble::iterator it;
-		}
+        // for (it = this->totals.begin(); it != this->totals.end(); it++) {
-	}
+        //   if (it->second < 1e-25) {
+        //     it->second = 0.0;
+        //   }
+        // }
 }
 void
 cxxSolution::Update_activities(const cxxNameDouble &original_tot)

src/phreeqcpp/Solution.h

@@ -105,13 +105,6 @@ class cxxSolution:public cxxNumKeyword
 	void dump_xml(std::ostream & os, unsigned int indent = 0) const;
 	void dump_raw(std::ostream & s_oss, unsigned int indent, int *n_out=NULL) const;
 
-        void dump_essential_names(std::vector<std::string> &e_names) const;
-
-        void get_essential_values(std::vector<LDBLE> &e_values,
-                                  const std::vector<std::string> &order) const;
-        void set_essential_values(std::vector<LDBLE>::iterator &it,
-                                  const std::vector<std::string> &order);
-
         LDBLE Get_master_activity(char *string) const;
 	void Set_master_activity(char *string, LDBLE value);	
 	LDBLE Get_total(const char *string) const;

src/phreeqcpp/Surface.cxx

@@ -178,105 +178,6 @@ cxxSurface::dump_raw(std::ostream & s_oss, unsigned int indent, int *n_out) cons
 	return;
 }
 
-void cxxSurface::dump_essential_names(std::vector<std::string> &e_names) {
-  // this is just some reserve size, might be changed in the future
-  const int reserve_size =
-      surface_comps.size() * 4 + surface_charges.size() * 100;
-  e_names.clear();
-  e_names.reserve(reserve_size);
-
-  for (auto &comp : this->surface_comps) {
-    const std::string phase_name = comp.Get_formula();
-    e_names.push_back(phase_name + "_moles");
-    e_names.push_back(phase_name + "_la");
-    e_names.push_back(phase_name + "_cb");
-
-    for (const auto &tot : comp.Get_totals()) {
-      e_names.push_back(phase_name + "_" + tot.first);
-    }
-  }
-
-  for (auto &charge : this->surface_charges) {
-    const std::string component = charge.Get_name() + "_charge";
-    e_names.push_back(component + "_area");
-    e_names.push_back(component + "_grams");
-    e_names.push_back(component + "_cb");
-    e_names.push_back(component + "_mw");
-    e_names.push_back(component + "_la");
-
-    for (const auto &tot : charge.Get_diffuse_layer_totals()) {
-      e_names.push_back(component + "_tot_" + tot.first);
-    }
-    // for (const auto &dl_layer : charge.Get_dl_species_map()) {
-    // 	e_names.push_back(component + "_dls_" +
-    // std::to_string(dl_layer.first));
-    // }
-  }
-}
-
-void cxxSurface::get_essential_values(std::vector<LDBLE> &e_values) {
-  // this is just some reserve size, might be changed in the future
-  const int reserve_size =
-      surface_comps.size() * 4 + surface_charges.size() * 100;
-  e_values.clear();
-  e_values.reserve(reserve_size);
-
-  for (auto &comp : this->surface_comps) {
-    e_values.push_back(comp.Get_moles());
-    e_values.push_back(comp.Get_la());
-    e_values.push_back(comp.Get_charge_balance());
-
-    for (const auto &tot : comp.Get_totals()) {
-      e_values.push_back(tot.second);
-    }
-  }
-
-  for (auto &charge : this->surface_charges) {
-    e_values.push_back(charge.Get_specific_area());
-    e_values.push_back(charge.Get_grams());
-    e_values.push_back(charge.Get_charge_balance());
-    e_values.push_back(charge.Get_mass_water());
-    e_values.push_back(charge.Get_la_psi());
-
-    for (const auto &tot : charge.Get_diffuse_layer_totals()) {
-      e_values.push_back(tot.second);
-    }
-    // for (const auto &dl_layer : charge.Get_dl_species_map()) {
-    // 	e_values.push_back(dl_layer.second);
-    // }
-  }
-}
-
-void cxxSurface::set_essential_values(std::vector<LDBLE>::iterator &it) {
-  for (auto &comp : this->surface_comps) {
-    comp.Set_moles(*(it++));
-    comp.Set_la(*(it++));
-    comp.Set_charge_balance(*(it++));
-
-    for (auto &tot : comp.Get_totals()) {
-      tot.second = *(it++);
-    }
-  }
-
-  for (auto &charge : this->surface_charges) {
-    charge.Set_specific_area(*(it++));
-    charge.Set_grams(*(it++));
-    charge.Set_charge_balance(*(it++));
-    charge.Set_mass_water(*(it++));
-    charge.Set_la_psi(*(it++));
-
-    for (auto &tot : charge.Get_diffuse_layer_totals()) {
-      tot.second = *(it++);
-    }
-
-    charge.Get_dl_species_map().clear();
-
-    // for (auto &dl_layer : charge.Get_dl_species_map()) {
-    // 	dl_layer.second = *(it++);
-    // }
-  }
-}
-
 void
 cxxSurface::read_raw(CParser & parser, bool check)
 {

src/phreeqcpp/Surface.h

@@ -1,3 +1,4 @@
+#include <set>
 #if !defined(SURFACE_H_INCLUDED)
 #define SURFACE_H_INCLUDED
 
@@ -31,11 +32,6 @@ public:
 	//void dump_xml(std::ostream & os, unsigned int indent = 0) const;
 	void dump_raw(std::ostream & s_oss, unsigned int indent, int *n_out=NULL) const;
 
-        void dump_essential_names(std::vector<std::string> &e_names);
-
-        void get_essential_values(std::vector<LDBLE> &e_values);
-        void set_essential_values(std::vector<LDBLE>::iterator &it);
-
         void read_raw(CParser & parser, bool check = true);
 	bool Get_related_phases(void) const;
 	bool Get_related_rate(void) const;

src/phreeqcpp/cxxKinetics.cxx

@@ -159,45 +159,6 @@ cxxKinetics::dump_raw(std::ostream & s_oss, unsigned int indent, int * n_out) co
 	return;
 }
 
-void cxxKinetics::dump_essential_names(
-    std::vector<std::string> &e_names) const {
-  e_names.clear();
-  e_names.reserve(this->kinetics_comps.size() * 2);
-
-  for (const auto &comp : this->kinetics_comps) {
-    const std::string name = comp.Get_rate_name();
-    e_names.push_back(name);
-    for (int i = 0; i < comp.Get_d_params().size(); i++) {
-      e_names.push_back(name + "_p" + std::to_string(i + 1));
-    }
-  }
-}
-
-void cxxKinetics::get_essential_values(std::vector<LDBLE> &e_values) const {
-  e_values.clear();
-  e_values.reserve(this->kinetics_comps.size() * 2);
-
-  for (const auto &comp : this->kinetics_comps) {
-    e_values.push_back(comp.Get_m());
-    for (const auto &param : comp.Get_d_params()) {
-      e_values.push_back(param);
-    }
-  }
-
-  return;
-}
-
-void cxxKinetics::set_essential_values(std::vector<LDBLE>::iterator &it) {
-  for (auto &comp : this->kinetics_comps) {
-    comp.Set_m(*(it++));
-    for (auto &param : comp.Get_d_params()) {
-      param = *(it++);
-    }
-  }
-
-  return;
-}
-
 void
 cxxKinetics::read_raw(CParser & parser, bool check)
 {

src/phreeqcpp/cxxKinetics.h

@@ -25,11 +25,6 @@ class cxxKinetics:public cxxNumKeyword
 
 	void dump_raw(std::ostream & s_oss, unsigned int indent, int * n_out=NULL) const;
 
-        void dump_essential_names(std::vector<std::string> &e_names) const;
-
-        void get_essential_values(std::vector<LDBLE> &e_values) const;
-        void set_essential_values(std::vector<LDBLE>::iterator &it);
-
         void read_raw(CParser & parser, bool check = true);
 
 	std::vector < LDBLE > &Get_steps(void) {return steps;}

src/phreeqcpp/mainsubs.cpp

@@ -341,6 +341,107 @@ set_use(void)
 	*/
 	return (OK);
 }
+
+int Phreeqc::initial_solutions_poet(int sol_id)
+/* ---------------------------------------------------------------------- */
+{
+  /*
+   *   Go through list of solutions, make initial solution calculations
+   *   for any marked "new".
+   */
+  int converge, converge1;
+  int last, n_user, print1;
+  char token[2 * MAX_LENGTH];
+
+  //   state = INITIAL_SOLUTION;
+  set_use();
+  print1 = TRUE;
+  dl_type_x = cxxSurface::NO_DL;
+  // std::map<int, cxxSolution>::iterator it = Rxn_solution_map.begin();
+  // for ( ; it != Rxn_solution_map.end(); it++)
+  //{
+  // for (size_t nn = 0; nn < Rxn_new_solution.size(); nn++)
+  cxxSolution &solution_ref = Rxn_solution_map.find(sol_id)->second;
+  initial_solution_isotopes = FALSE;
+  //   if (solution_ref.Get_new_def()) {
+  use.Set_solution_ptr(&solution_ref);
+  LDBLE d0 = solution_ref.Get_density();
+  // LDBLE d1 = 0;
+  bool diag = (diagonal_scale == TRUE) ? true : false;
+  int count_iterations = 0;
+  //   std::string input_units = solution_ref.Get_initial_data()->Get_units();
+  //   cxxISolution *initial_data_ptr = solution_ref.Get_initial_data();
+  density_iterations = 0;
+  // for (;;) {
+  prep();
+  k_temp(solution_ref.Get_tc(), solution_ref.Get_patm());
+  set(TRUE);
+  always_full_pitzer = FALSE;
+
+  diagonal_scale = TRUE;
+  converge = model();
+  if (converge == ERROR /*&& diagonal_scale == FALSE*/) {
+    diagonal_scale = TRUE;
+    always_full_pitzer = TRUE;
+    set(TRUE);
+    converge = model();
+  }
+  //   density_iterations++;
+  //   if (solution_ref.Get_initial_data()->Get_calc_density()) {
+  //     solution_ref.Set_density(calc_dens());
+  //     if (!equal(d0, solution_ref.Get_density(), 1e-8)) {
+  //       initial_data_ptr->Set_units(input_units);
+  //       d0 = solution_ref.Get_density();
+  //       if (count_iterations++ < 20) {
+  //         diag = (diagonal_scale == TRUE) ? true : false;
+  //         continue;
+  //       } else {
+  //         error_msg(
+  //             sformatf("%s %d.",
+  //                      "Density calculation failed for initial solution ",
+  //                      solution_ref.Get_n_user()),
+  //             STOP);
+  //       }
+  //     }
+  //   }
+  //   break;
+  // }
+  diagonal_scale = (diag) ? TRUE : FALSE;
+  converge1 = check_residuals();
+  sum_species();
+  viscosity();
+  add_isotopes(solution_ref);
+  punch_all();
+  //   print_all();
+  density_iterations = 0;
+  /* free_model_allocs(); */
+  // remove pr_in
+  for (size_t i = 0; i < count_unknowns; i++) {
+    if (x[i]->type == SOLUTION_PHASE_BOUNDARY)
+      x[i]->phase->pr_in = false;
+  }
+
+  if (converge == ERROR || converge1 == ERROR) {
+    error_msg(sformatf("%s %d.",
+                       "Model failed to converge for initial solution ",
+                       solution_ref.Get_n_user()),
+              STOP);
+  }
+  n_user = solution_ref.Get_n_user();
+  last = solution_ref.Get_n_user_end();
+  /* copy isotope data */
+  if (solution_ref.Get_isotopes().size() > 0) {
+    isotopes_x = solution_ref.Get_isotopes();
+  } else {
+    isotopes_x.clear();
+  }
+  xsolution_save(n_user);
+  Utilities::Rxn_copies(Rxn_solution_map, n_user, last);
+  //   }
+  initial_solution_isotopes = FALSE;
+  return (OK);
+}
+
 /* ---------------------------------------------------------------------- */
 int Phreeqc::
 initial_solutions(int print)

src/phreeqcpp/model.cpp

@@ -5003,16 +5003,16 @@ surface_model(void)
  */
 	debug_diffuse_layer_save = debug_diffuse_layer;
 	debug_model_save = debug_model;
-	if (last_model.force_prep)
+        // if (last_model.force_prep)
-	{
+        // {
+        // 	same_model = FALSE;
+        // }
+        // else
+        // {
+        // 	same_model = check_same_model();
+        // }
 		same_model = FALSE;
-	}
+        if (dl_type_x != cxxSurface::NO_DL && same_model == FALSE) {
-	else
-	{
-		same_model = check_same_model();
-	}
-	if (dl_type_x != cxxSurface::NO_DL && same_model == FALSE)
-	{
 		s_diff_layer.clear();
 		for (i = 0; i < (int)s.size(); i++)
 		{

src/phreeqcpp/prep.cpp

@@ -36,19 +36,9 @@ prep(void)
  */
 	cxxSolution *solution_ptr;
 
-	if (state >= REACTION)
+        // UP: we force to reset the model
-	{
-		same_model = check_same_model();
-	}
-	else
-	{
         same_model = FALSE;
         last_model.force_prep = true;
-	}
-	/*same_model = FALSE; */
-/*
- *   Initialize s, master, and unknown pointers
- */
         solution_ptr = use.Get_solution_ptr();
         if (solution_ptr == NULL)
 	{
@@ -87,10 +77,15 @@ prep(void)
  *   Allocate space for array
  */
 		my_array.resize((max_unknowns + 1) * max_unknowns);
+                for (double &val : my_array) {
+                  val = 0;
+                }
                 delta.resize(max_unknowns);
+                for (double &val : delta) {
+                  val = 0;
+                }
                 residual.resize(max_unknowns);
-		for (int j = 0; j < max_unknowns; j++)
+                for (int j = 0; j < max_unknowns; j++) {
-		{
                   residual[j] = 0;
                 }