feat: add PhreeqcRunner class and integrate with PhreeqcMatrix; refactor constructors and update CMakeLists.txt

poet/CMakeLists.txt

@@ -1,5 +1,6 @@
 set(POET_SOURCE_FILES
-    src/Engine.cpp 
+    src/Engine.cpp
+    src/Runner.cpp 
     #Wrappers
     src/Wrapper/EquilibriumWrapper.cpp
     src/Wrapper/EquilibriumCompWrapper.cpp
@@ -33,6 +34,7 @@ if (BUILD_TESTING AND STANDALONE_BUILD)
     set(POET_TEST_SOURCE_FILES
         test/testPhreeqcEngine.cpp
         test/testPhreeqcMatrix.cpp
+        test/testPhreeqcRunner.cpp
         test/utils.cpp
     )
 

poet/include/PhreeqcMatrix.hpp

@@ -57,7 +57,7 @@ public:
    * object!
    * @param other
    */
-  PhreeqcMatrix(const PhreeqcMatrix &other);
+  PhreeqcMatrix(const PhreeqcMatrix &other) = default;
 
   /**
    * @brief Construct a new Phreeqc Matrix object
@@ -66,7 +66,7 @@ public:
    * object!
    * @param other
    */
-  PhreeqcMatrix(PhreeqcMatrix &&other);
+  PhreeqcMatrix(PhreeqcMatrix &&other) = default;
 
   /**
    * @brief Assignment operator
@@ -75,7 +75,7 @@ public:
    * @param other
    * @return PhreeqcMatrix&
    */
-  PhreeqcMatrix &operator=(const PhreeqcMatrix &other);
+  PhreeqcMatrix &operator=(const PhreeqcMatrix &other) = default;
 
   /**
    * @brief Assignment operator
@@ -84,7 +84,7 @@ public:
    * @param other
    * @return PhreeqcMatrix&
    */
-  PhreeqcMatrix &operator=(PhreeqcMatrix &&other);
+  PhreeqcMatrix &operator=(PhreeqcMatrix &&other) = default;
 
   /**
    * @brief Access the value of a given cell by name.

poet/include/PhreeqcRunner.hpp

@@ -0,0 +1,66 @@
+#pragma once
+
+#include "PhreeqcEngine.hpp"
+#include "PhreeqcMatrix.hpp"
+#include <cstddef>
+#include <memory>
+#include <unordered_map>
+#include <vector>
+
+/**
+ * @class PhreeqcRunner
+ * @brief Manages the execution of Phreeqc simulations.
+ *
+ * The PhreeqcRunner class is responsible for managing and running Phreeqc
+ * simulations. It initializes a set of PhreeqcEngine instances based on the
+ * provided PhreeqcMatrix and provides functionality to run simulations and
+ * retrieve the number of engines.
+ *
+ * @note Copy and move operations are deleted to prevent unintended behavior.
+ */
+class PhreeqcRunner {
+public:
+  PhreeqcRunner(const PhreeqcRunner &) = delete;
+  PhreeqcRunner(PhreeqcRunner &&) = delete;
+  PhreeqcRunner &operator=(const PhreeqcRunner &) = delete;
+  PhreeqcRunner &operator=(PhreeqcRunner &&) = delete;
+
+  /**
+   * @brief Constructs a PhreeqcRunner object with the given PhreeqcMatrix.
+   *
+   * Constructs a PhreeqcRunner object with the given PhreeqcMatrix. For each
+   *  cell found in the PhreeqcMatrix, a PhreeqcEngine is created and stored in
+   *  an internal map.
+   *
+   * @param matrix A reference to a PhreeqcMatrix object used to initialize the
+   * PhreeqcRunner.
+   */
+  PhreeqcRunner(const PhreeqcMatrix &matrix);
+  ~PhreeqcRunner() = default;
+
+  /**
+   * @brief Runs the simulation with the given input and output data for a
+   * specified time step.
+   *
+   * @param simulationInOut A reference to a 2D vector containing the input data
+   * for the simulation. The vector will be modified to contain the output data
+   * after the simulation.
+   * @param time_step The time step for the simulation.
+   */
+  void run(std::vector<std::vector<double>> &simulationInOut,
+           const double time_step);
+
+  /**
+   * @brief Returns the number of engines currently stored.
+   *
+   * This function provides the count of engines that are currently
+   * stored in the _engineStorage container.
+   *
+   * @return std::size_t The number of engines.
+   */
+  std::size_t numEngines() const { return _engineStorage.size(); }
+
+private:
+  std::unordered_map<int, std::unique_ptr<PhreeqcEngine>> _engineStorage;
+  std::vector<double> _buffer;
+};

poet/src/Engine.cpp

@@ -110,8 +110,8 @@ void PhreeqcEngine::runCell(std::vector<double> &cell_values,
     throw std::invalid_argument("Time step must be positive");
   }
 
-  // skip ID
-  std::span<double> cell_data{cell_values.begin() + 1, cell_values.end()};
+  // ID is already skipped by PhreeqcRunner, so no need to start ahead
+  std::span<double> cell_data{cell_values.begin(), cell_values.end()};
 
   this->impl->set_essential_values(cell_data);
   this->impl->run(time_step);

poet/src/PhreeqcMatrix/Access.cpp

@@ -7,9 +7,10 @@
 #include <vector>
 
 PhreeqcMatrix PhreeqcMatrix::subset(const std::vector<int> &indices) const {
-  PhreeqcMatrix result;
+  PhreeqcMatrix result(*this);
 
-  result._m_pqc = _m_pqc;
+  result._m_map.clear();
+  result._m_internal_names.clear();
 
   for (const auto &index : indices) {
     result._m_map[index] = _m_map.at(index);

poet/src/PhreeqcMatrix/Ctor.cpp

@@ -17,32 +17,32 @@ PhreeqcMatrix::PhreeqcMatrix(const std::string &database,
   this->initialize();
 }
 
-PhreeqcMatrix::PhreeqcMatrix(const PhreeqcMatrix &other)
-    : _m_map(other._m_map), _m_internal_names(other._m_internal_names),
-      _m_surface_primaries(other._m_surface_primaries), _m_pqc(other._m_pqc),
-      _m_database(other._m_database) {}
+// PhreeqcMatrix::PhreeqcMatrix(const PhreeqcMatrix &other)
+//     : _m_map(other._m_map), _m_internal_names(other._m_internal_names),
+//       _m_surface_primaries(other._m_surface_primaries), _m_pqc(other._m_pqc),
+//       _m_database(other._m_database) {}
 
-PhreeqcMatrix::PhreeqcMatrix(PhreeqcMatrix &&other)
-    : _m_map(other._m_map), _m_internal_names(other._m_internal_names),
-      _m_surface_primaries(other._m_surface_primaries), _m_pqc(other._m_pqc),
-      _m_database(other._m_database) {}
+// PhreeqcMatrix::PhreeqcMatrix(PhreeqcMatrix &&other)
+//     : _m_map(other._m_map), _m_internal_names(other._m_internal_names),
+//       _m_surface_primaries(other._m_surface_primaries), _m_pqc(other._m_pqc),
+//       _m_database(other._m_database) {}
 
-PhreeqcMatrix &PhreeqcMatrix::operator=(const PhreeqcMatrix &other) {
-  _m_map = other._m_map;
-  _m_internal_names = other._m_internal_names;
-  _m_surface_primaries = other._m_surface_primaries;
-  _m_pqc = other._m_pqc;
-  _m_database = other._m_database;
+// PhreeqcMatrix &PhreeqcMatrix::operator=(const PhreeqcMatrix &other) {
+//   _m_map = other._m_map;
+//   _m_internal_names = other._m_internal_names;
+//   _m_surface_primaries = other._m_surface_primaries;
+//   _m_pqc = other._m_pqc;
+//   _m_database = other._m_database;
 
-  return *this;
-}
+//   return *this;
+// }
 
-PhreeqcMatrix &PhreeqcMatrix::operator=(PhreeqcMatrix &&other) {
-  _m_map = other._m_map;
-  _m_internal_names = other._m_internal_names;
-  _m_surface_primaries = other._m_surface_primaries;
-  _m_pqc = other._m_pqc;
-  _m_database = other._m_database;
+// PhreeqcMatrix &PhreeqcMatrix::operator=(PhreeqcMatrix &&other) {
+//   _m_map = other._m_map;
+//   _m_internal_names = other._m_internal_names;
+//   _m_surface_primaries = other._m_surface_primaries;
+//   _m_pqc = other._m_pqc;
+//   _m_database = other._m_database;
 
-  return *this;
-}
+//   return *this;
+// }

poet/src/Runner.cpp

@@ -0,0 +1,57 @@
+#include "PhreeqcEngine.hpp"
+#include "PhreeqcMatrix.hpp"
+#include "PhreeqcRunner.hpp"
+#include <cmath>
+#include <cstddef>
+#include <memory>
+#include <vector>
+
+PhreeqcRunner::PhreeqcRunner(const PhreeqcMatrix &matrix) {
+  // first make sure to have enough space in our buffer
+  this->_buffer.reserve(matrix.get().names.size());
+
+  // Create a PhreeqcEngine for each id
+  for (const auto &id : matrix.getIds()) {
+    this->_engineStorage[id] = std::make_unique<PhreeqcEngine>(matrix, id);
+  }
+}
+
+static void copy_to_buffer(std::vector<double> &buffer,
+                           const std::vector<double> &current_inout) {
+  for (std::size_t j = 1; j < current_inout.size(); j++) {
+    if (std::isnan(current_inout[j])) {
+      continue;
+    }
+    buffer.push_back(current_inout[j]);
+  }
+}
+
+static void copy_from_buffer(const std::vector<double> &buffer,
+                             std::vector<double> &current_inout) {
+  std::size_t buffer_index = 0;
+  for (std::size_t j = 1; j < current_inout.size(); j++) {
+    if (std::isnan(current_inout[j])) {
+      continue;
+    }
+    current_inout[j] = buffer[buffer_index++];
+  }
+}
+
+void PhreeqcRunner::run(std::vector<std::vector<double>> &simulationInOut,
+                        const double time_step) {
+  for (std::size_t i = 0; i < simulationInOut.size(); i++) {
+    const auto &current_inout = simulationInOut[i];
+    const auto pqc_id = static_cast<int>(current_inout[0]);
+
+    this->_buffer.clear();
+
+    // Copy the input to the buffer while ignoring the first element and NaNs
+    copy_to_buffer(this->_buffer, current_inout);
+
+    this->_engineStorage.at(pqc_id)->runCell(this->_buffer, time_step);
+
+    // Copy the buffer back to the output while ignoring the first element and
+    // NaNs
+    copy_from_buffer(this->_buffer, simulationInOut[i]);
+  }
+}

poet/test/testInput.hpp.in

@@ -4,6 +4,8 @@
 #include <string>
 #include <vector>
 
+#include <gtest/gtest.h>
+
 #define POET_TEST(name) TEST(TestPOET, name)
 
 namespace base_test {

poet/test/testPhreeqcRunner.cpp

@@ -0,0 +1,75 @@
+#include "PhreeqcRunner.hpp"
+#include "utils.hpp"
+
+#include <cmath>
+#include <cstddef>
+#include <gtest/gtest.h>
+#include <testInput.hpp>
+#include <vector>
+
+const std::string test_database = readFile(barite_test::database);
+const std::string test_script = readFile(barite_test::script);
+
+constexpr std::size_t num_cells = 10;
+
+POET_TEST(PhreeqcRunnerConstructor) {
+  PhreeqcMatrix pqc_mat(test_database, test_script);
+  EXPECT_NO_THROW(PhreeqcRunner tmp(pqc_mat));
+}
+
+POET_TEST(PhreeqcRunnerSimulation) {
+  PhreeqcMatrix pqc_mat(test_database, test_script);
+  const auto subsetted_pqc_mat = pqc_mat.subset({2, 3});
+  PhreeqcRunner runner(subsetted_pqc_mat);
+
+  const auto stl_mat = subsetted_pqc_mat.get();
+  const auto matrix_values = stl_mat.values;
+  const auto num_columns = stl_mat.names.size();
+
+  std::vector<std::vector<double>> simulationInOut;
+
+  for (std::size_t index = 0; index < num_cells; ++index) {
+    if (index < static_cast<std::size_t>(num_cells / 2)) {
+      simulationInOut.push_back(std::vector<double>(
+          matrix_values.begin(), matrix_values.begin() + num_columns));
+    } else {
+      simulationInOut.push_back(std::vector<double>(
+          matrix_values.begin() + num_columns, matrix_values.end()));
+    }
+  }
+  EXPECT_NO_THROW(runner.run(simulationInOut, 100));
+
+  constexpr std::size_t half_cells = num_cells / 2;
+  constexpr int expected_value_first_half = 2;
+  constexpr int expected_value_second_half = 3;
+
+  for (std::size_t cell_index = 0; cell_index < simulationInOut.size();
+       ++cell_index) {
+    const bool is_first_half = cell_index < half_cells;
+    if (is_first_half) {
+      EXPECT_EQ(simulationInOut[cell_index][0], expected_value_first_half);
+      EXPECT_TRUE(std::isnan(simulationInOut[cell_index][11]));
+    } else {
+      EXPECT_EQ(simulationInOut[cell_index][0], expected_value_second_half);
+      EXPECT_FALSE(std::isnan(simulationInOut[cell_index][11]));
+    }
+
+    EXPECT_NEAR(simulationInOut[cell_index][1], 111, 1);
+    EXPECT_NEAR(simulationInOut[cell_index][2], 55, 1);
+  }
+}
+
+POET_TEST(PhreeqcRunnerUnknownID) {
+  PhreeqcMatrix pqc_mat(test_database, test_script);
+  PhreeqcRunner runner(pqc_mat);
+
+  std::vector<std::vector<double>> simulationInOut;
+
+  for (std::size_t index = 0; index < num_cells; ++index) {
+    simulationInOut.push_back(std::vector<double>(pqc_mat.get().names.size()));
+  }
+
+  simulationInOut[0][0] = 1000;
+
+  EXPECT_THROW(runner.run(simulationInOut, 100), std::out_of_range);
+}