From f6f89a11614757683080e3ff9a12ad0c34b2184c Mon Sep 17 00:00:00 2001 From: Marco De Lucia Date: Mon, 2 Dec 2024 18:20:02 +0100 Subject: [PATCH] Fixed commit --- poet/test/barite_het.pqi | 13 ++- poet/test/testGolemRunner.cpp | 161 ++++++++++++++++++++++++++++++++++ 2 files changed, 167 insertions(+), 7 deletions(-) create mode 100644 poet/test/testGolemRunner.cpp diff --git a/poet/test/barite_het.pqi b/poet/test/barite_het.pqi index b8cf1e10..8d3b16e5 100644 --- a/poet/test/barite_het.pqi +++ b/poet/test/barite_het.pqi @@ -23,12 +23,11 @@ RUN_CELLS COPY solution 1 2-3 ## Here a 5x2 domain: - - |---+---+---+---+---| - -> | 2 | 2 | 2 | 2 | 2 | - 4 |---+---+---+---+---| - -> | 3 | 3 | 3 | 3 | 3 | - |---+---+---+---+---| +## |---+---+---+---+---| +## -> | 2 | 2 | 2 | 2 | 2 | +## 4 |---+---+---+---+---| +## -> | 3 | 3 | 3 | 3 | 3 | +## |---+---+---+---+---| ## East boundary: "injection" of solution 4. North, W, S boundaries: closed @@ -77,4 +76,4 @@ END ## KINETICS 3. RUN_CELLS - -cells 2-4 \ No newline at end of file + -cells 2-4 diff --git a/poet/test/testGolemRunner.cpp b/poet/test/testGolemRunner.cpp new file mode 100644 index 00000000..25553ace --- /dev/null +++ b/poet/test/testGolemRunner.cpp @@ -0,0 +1,161 @@ +// Time-stamp: "Last modified 2024-12-02 17:37:08 delucia" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "PhreeqcMatrix.hpp" +#include "PhreeqcEngine.hpp" +#include "PhreeqcRunner.hpp" + + +std::string readFile(const std::string &path) { + std::string string_rpath(PATH_MAX, '\0'); + + if (realpath(path.c_str(), string_rpath.data()) == nullptr) { + throw std::runtime_error(":: Failed to resolve the realpath to file " + path); + } + + std::ifstream file(string_rpath); + + if (!file.is_open()) { + throw std::runtime_error(":: Failed to open file: " + path); + } + + std::stringstream buffer; + buffer << file.rdbuf(); + + return buffer.str(); +} + +// pretty print a vector, standard implementation from stackoverflow +template +std::ostream & operator<<(std::ostream & os, std::vector vec) +{ + os << "{ "; + std::copy(vec.begin(), vec.end(), std::ostream_iterator(os, ", ")); + os << " }"; + return os; +} + + +int main(int argc, char *argv[]) { + + if (argc < 3) { + std::cout << "::" << argv[0] << ": two args needed, script and database\n"; + return 1; + } + + + ////// INITIALISATION + // read Script and Database and put it in a std::string + auto script = readFile(argv[1]); + auto db = readFile(argv[2]); + + // Create the matrix directly from database and init script + PhreeqcMatrix pqc_mat(db, script); + + // How many different SOLUTIONS ("CELLS") are defined in the script? + const auto ids = pqc_mat.getIds(); + + int n = ids.size(); + + std::cout << ":: Found " << n << " distinct PHREEQC problems \n"; + std::cout << ids << "\n"; + + const auto solutes = pqc_mat.getSolutionNames(); + std::cout << ":: These are the common solutes across all the " << n << " problems: \n"; + std::cout << solutes << "\n"; + + // iterate on the ids (THEY start at 1!!) + for (const auto & i : ids) { + auto pphases = pqc_mat.getEquilibriumNames(i); + if (!pphases.empty()) { + std::cout << ":: Equilibrium phases [" << (int) i << "]: \n"; + std::cout << pphases << "\n"; + } + + auto kinetics = pqc_mat.getKineticsNames(i); + if (!kinetics.empty()) { + std::cout << ":: Kinetics [" << i << "]: \n"; + std::cout << kinetics << "\n"; + } + } + + // The exported data type holds the matrix in a "STL format" with + // a "header" of names and their accompanying values. The values + // are stored in a row-major order per default. + auto exported_mat = pqc_mat.get(); + // Get the total number of solutes + const int len = exported_mat.names.size(); + // Get the values as reference to modify them in place + std::vector &cell_values = exported_mat.values; + + std::cout << ":: Values in the PhreeqcMatrix: \n"; + + // std::cout << exported_mat.names << "\n"; + // std::cout << cell_values << "\n"; + // END INIT + + + + //// Phreeqc RUN through the new Runner class + + // optional SUBSET the matrix (i.e., the unique ids defined in + // golem map as input) + // const auto subsetted_pqc_mat = pqc_mat.subset({1, 2}); + PhreeqcRunner runner(pqc_mat); + + const auto stl_mat = pqc_mat.get(); + const auto matrix_values = stl_mat.values; + const auto num_columns = stl_mat.names.size(); + const auto spec_names = stl_mat.names; + + // container to pass in/out + std::vector> simulationInOut; + + // grid cells + const std::size_t num_cells = 10; + const std::size_t half_cells = 5; + + // copy the values to the InOut vector. We replicate cell 1 + for (std::size_t index = 0; index < num_cells; ++index) { + if (index < half_cells) { + simulationInOut.push_back(std::vector( + matrix_values.begin(), matrix_values.begin() + num_columns)); + } else { + simulationInOut.push_back(std::vector( + matrix_values.begin() + num_columns, matrix_values.end())); + } + } + + const double timestep = 100.; + + // compute 1 timestep + runner.run(simulationInOut, timestep); + + for (std::size_t cell_index = 0; cell_index < simulationInOut.size(); ++cell_index) { + const bool is_first_half = cell_index < half_cells; + if (is_first_half) { + std::cout << "Grid element: " << cell_index << " \n"; + for (std::size_t spec = 0; spec < num_columns; ++spec) { + std::cout << ":" << spec_names[spec] << "=" << simulationInOut[cell_index][spec]; + } + std::cout << "\n"; + } + } + + + return 0; +} + +// Oneliner for rz-vm278 relative to iphreeqc/poet/test!! + +// g++ testGolemRunner.cpp -o testG -Wall -I../../poet/include -I../../src -I../../src/phreeqcpp -I../../src/phreeqcpp/common -I../../src/phreeqcpp/PhreeqcKeywords -lIPhreeqc -lIPhreeqcPOET -L../../bbuild/ -L../../bbuild/poet