Merge branch 'ml/remove-cpp20-api' into 'poet'

Improve interoperability with exclusive Cpp language features

See merge request naaice/iphreeqc!9
This commit is contained in:
Max Lübke 2024-09-10 18:18:56 +02:00
commit e02969c891
27 changed files with 2242 additions and 189 deletions

View File

@ -393,6 +393,7 @@ if (BUILD_CLR_LIBS)
endif() endif()
add_subdirectory(poet) add_subdirectory(poet)
add_subdirectory(app)
include (CTest) include (CTest)

25
app/CMakeLists.txt Normal file
View File

@ -0,0 +1,25 @@
## Time-stamp: "Last modified 2024-08-26 17:42:18 delucia"
add_executable(ipqMain
ipqMain.cpp
)
add_executable(ipqEngine
ipqEngine.cpp
)
# add_library(IPhreeqcPOET ${PQC_POET_SRC})
##target_link_libraries(IPhreeqcPOET PUBLIC IPhreeqc)
target_include_directories(IPhreeqcPOET PUBLIC
$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
$<INSTALL_INTERFACE:include>
PRIVATE src)
target_link_libraries(ipqMain
IPhreeqcPOET
)
target_link_libraries(ipqEngine
IPhreeqcPOET
)

140
app/ipqEngine.cpp Normal file
View File

@ -0,0 +1,140 @@
#include <algorithm>
#include <cstddef>
#include <cstdio>
#include <iomanip>
#include <linux/limits.h>
#include <stdlib.h>
#include <string>
#include <unistd.h>
#include <vector>
#include <wordexp.h>
#include "POETInit.hpp"
#include "PhreeqcEngine.hpp"
#include "PhreeqcInit.hpp"
static 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();
}
void printVectorsAsTable(const std::vector<std::string> &names,
const std::vector<double> &values) {
if (names.size() != values.size()) {
std::cerr << "Error: Vectors are not of the same length." << std::endl;
return;
}
std::cout << std::left << std::setw(20) << "Array Index" << std::setw(20)
<< "Name" << std::setw(20) << "Value" << std::endl;
std::cout << std::string(40, '-') << std::endl;
for (std::size_t i = 0; i < names.size(); ++i) {
std::cout << std::left << std::setw(20) << std::to_string(i)
<< std::setw(20) << names[i] << std::setw(20) << values[i]
<< std::endl;
}
}
template <class T> void printVector(const std::vector<T> &vec) {
for (const auto &elem : vec) {
std::cout << elem << " ";
}
std::cout << std::endl;
}
static std::string TestInput = R"(SOLUTION 1
units mol/kgw
temp 25
Ca 0.1
Mg 0.1
Cl 0.5 charge
Na 0.1
PURE 1
Calcite 0.0 1
Dolomite 0.0 0
## RUN_CELLS
## -cells 1
END)";
int main(int argc, char *argv[]) {
if (argc < 2) {
std::cout << "PHREEQC database needed as argument" << std::endl;
return 1;
}
const std::string &database = argv[1];
wordexp_t exp_result;
// std::cout << "Reading script from file " << script << " using database "
// << database << std::endl;
// expand the script and database path
// const std::string read_script = readFile(script);
const std::string read_database = readFile(database);
// initialize the module instance "ipqcmod"
PhreeqcInit ipqcmod(read_database, TestInput);
// allocate tabular data
PhreeqcInit::PhreeqcMat pqc_mat = ipqcmod.getPhreeqcMat();
// for (std::size_t i = 0; i < nrows; ++i) {
auto pure_names = ipqcmod.getEquilibriumNames(pqc_mat.ids[0]);
auto totals_names = ipqcmod.getSolutionNames(pqc_mat.ids[0]);
POETConfig MyConfig;
MyConfig.database = read_database;
MyConfig.input_script = TestInput;
MyConfig.cell = POETInitCell{
// H, O and charge need to be removed. I'm not sure why I did it this way.
std::vector<std::string>(totals_names.begin() + 3, totals_names.end()),
{},
{},
{},
pure_names,
{},
{}};
std::vector<std::string> full_names = {"ID"};
full_names.insert(full_names.end(), pqc_mat.names.begin(),
pqc_mat.names.end());
PhreeqcEngine MyEngine(MyConfig);
std::vector<double> NewInp = pqc_mat.values[0];
NewInp.insert(NewInp.begin(), 1);
std::cout << "Before running the cell:\n";
printVectorsAsTable(full_names, NewInp);
// std::for_each(NewInp.begin()+3, NewInp.end(),
// [](double a) { return a * 0.8;});
NewInp[5] += 0.2; // Cl
NewInp[7] += 0.2; // Na
MyEngine.runCell(NewInp, 200);
std::cout << "\n\nAfter running the cell:\n";
printVectorsAsTable(full_names, NewInp);
return 0;
}

92
app/ipqMain.cpp Normal file
View File

@ -0,0 +1,92 @@
// Time-stamp: "Last modified 2024-08-27 09:00:18 delucia"
#include <cstddef>
#include <linux/limits.h>
#include <stdlib.h>
#include <string>
#include <unistd.h>
#include <vector>
#include <wordexp.h>
#include "POETInit.hpp"
#include "PhreeqcEngine.hpp"
#include "PhreeqcInit.hpp"
static 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();
}
int main(int argc, char *argv[]) {
const std::string &script = argv[1];
const std::string &database = argv[2];
// wordexp_t exp_result;
std::cout << "Reading script from file " << script << " using database "
<< database << std::endl;
// expand the script and database path
const std::string read_script = readFile(script);
const std::string read_database = readFile(database);
// initialize the module instance "ipqcmod"
PhreeqcInit ipqcmod(read_database, read_script);
// allocate tabular data
PhreeqcInit::PhreeqcMat pqc_mat = ipqcmod.getPhreeqcMat();
const std::size_t ncols = pqc_mat.names.size();
const std::size_t nrows = pqc_mat.values.size();
std::cout << "Script contains " << nrows << " rows and " << ncols
<< " columns" << std::endl;
// pqc_mat.names.insert(pqc_mat.names.begin(), "ID");
// Rcpp::NumericMatrix mat(nrows, ncols + 1);
std::cout << std::endl << "Number of CELLS/SOLUTIONS: " << nrows << std::endl;
for (std::size_t i = 0; i < nrows; ++i) {
auto pure_names = ipqcmod.getEquilibriumNames(pqc_mat.ids[i]);
auto totals_names = ipqcmod.getSolutionNames(pqc_mat.ids[i]);
std::cout << std::endl << "SOL " << pqc_mat.ids[i] << " - Pure phases : ";
for (std::size_t j = 0; j < pure_names.size(); j++) {
std::cout << pure_names[j] << ", ";
}
std::cout << std::endl << "SOL " << pqc_mat.ids[i] << " - Total concs: ";
for (std::size_t j = 0; j < totals_names.size(); j++) {
std::cout << totals_names[j] << ", ";
}
std::cout << std::endl;
}
// alternatively, use the pqc_mat object
std::cout << std::endl << "Total pqc_mat: " << std::endl;
for (std::size_t i = 0; i < nrows; ++i) {
std::cout << "ID: " << pqc_mat.ids[i] << " - ";
for (std::size_t j = 0; j < ncols; ++j) {
std::cout << pqc_mat.names[j] << ": " << pqc_mat.values[i][j] << ", ";
}
std::cout << std::endl;
}
std::cout << std::endl;
return 0;
}

1853
app/phreeqc.dat Normal file

File diff suppressed because it is too large Load Diff

13
app/testEngine.pqi Normal file
View File

@ -0,0 +1,13 @@
SOLUTION 1
units mol/kgw
temp 25
Ca 0.1
Mg 0.1
Cl 0.5 charge
Na 0.1
PURE 1
Calcite 0.0 1
Dolomite 0.0 0
## RUN_CELLS
## -cells 1
END

View File

@ -9,4 +9,4 @@ target_include_directories(IPhreeqcPOET PUBLIC
PRIVATE src) PRIVATE src)
# Set C++23 standard # Set C++23 standard
target_compile_features(IPhreeqcPOET PUBLIC cxx_std_23) target_compile_features(IPhreeqcPOET PUBLIC cxx_std_20)

View File

@ -1,49 +1,35 @@
#pragma once #pragma once
#include "EquilibriumWrapper.hpp"
#include "ExchangeWrapper.hpp"
#include "KineticWrapper.hpp"
#include "POETInit.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 <memory> #include <memory>
#include <span>
#include <string>
#include <vector> #include <vector>
/** /**
* @brief Class for running Phreeqc wrappped in POET * @brief Class for running Phreeqc wrappped in POET
* *
* Direct interface to Phreeqc, without utilizing *_MODIFY keywords/scripts to * Direct interface to Phreeqc, without utilizing *_MODIFY keywords/scripts
to
* set new values. Use with already initialized Phreeqc config. * set new values. Use with already initialized Phreeqc config.
* *
*/ */
class PhreeqcEngine : public IPhreeqc { class PhreeqcEngine {
public: public:
/** /**
* @brief Construct a new Phreeqc Engine object * @brief Construct a new Phreeqc Engine object
* *
* Construct a new Phreeqc Engine object by previously initialized POETConfig. * Construct a new Phreeqc Engine object by previously initialized
POETConfig.
* *
* @param config Holds the configuration for the Phreeqc engine. * @param config Holds the configuration for the Phreeqc engine.
*/ */
PhreeqcEngine(const POETConfig &config) : IPhreeqc() { PhreeqcEngine(const POETConfig &config);
this->LoadDatabaseString(config.database.c_str());
this->RunString(config.input_script.c_str());
this->init_wrappers(config.cell); /**
} * @brief Destroy the Phreeqc Engine object
*
*/
~PhreeqcEngine();
/** /**
* @brief Siimulate a cell for a given time step * @brief Siimulate a cell for a given time step
@ -56,42 +42,6 @@ public:
void runCell(std::vector<double> &cell_values, double time_step); void runCell(std::vector<double> &cell_values, double time_step);
private: private:
void init_wrappers(const POETInitCell &cell); class Impl;
void run(double time_step); std::unique_ptr<Impl> impl;
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;
}; };

View File

@ -102,7 +102,7 @@ public:
* cell. Empty if no exchange is defined. * cell. Empty if no exchange is defined.
*/ */
std::vector<std::string> getExchanger(int id) { std::vector<std::string> getExchanger(int id) {
if (this->exchanger.contains(id)) { if (contains_id(this->exchanger, id)) {
return this->exchanger[id]; return this->exchanger[id];
} }
@ -117,7 +117,7 @@ public:
* cell. Empty if no kinetics are defined. * cell. Empty if no kinetics are defined.
*/ */
std::vector<std::string> getKineticsNames(int id) { std::vector<std::string> getKineticsNames(int id) {
if (this->kinetics.contains(id)) { if (contains_id(this->kinetics, id)) {
return this->kinetics[id]; return this->kinetics[id];
} }
@ -132,7 +132,7 @@ public:
* cell. Empty if no equilibrium is defined. * cell. Empty if no equilibrium is defined.
*/ */
std::vector<std::string> getEquilibriumNames(int id) { std::vector<std::string> getEquilibriumNames(int id) {
if (this->equilibrium.contains(id)) { if (contains_id(this->equilibrium, id)) {
return this->equilibrium[id]; return this->equilibrium[id];
} }
@ -147,7 +147,7 @@ public:
* for the cell. Empty if no surface is defined. * for the cell. Empty if no surface is defined.
*/ */
std::vector<std::string> getSurfaceCompNames(int id) { std::vector<std::string> getSurfaceCompNames(int id) {
if (this->surface_comps.contains(id)) { if (contains_id(this->surface_comps, id)) {
return this->surface_comps[id]; return this->surface_comps[id];
} }
@ -162,7 +162,7 @@ public:
* the cell. Empty if no surface is defined. * the cell. Empty if no surface is defined.
*/ */
std::vector<std::string> getSurfaceChargeNames(int id) { std::vector<std::string> getSurfaceChargeNames(int id) {
if (this->surface_charge.contains(id)) { if (contains_id(this->surface_charge, id)) {
return this->surface_charge[id]; return this->surface_charge[id];
} }
@ -216,6 +216,15 @@ private:
std::map<int, std::vector<std::string>> surface_comps; std::map<int, std::vector<std::string>> surface_comps;
std::map<int, std::vector<std::string>> surface_charge; std::map<int, std::vector<std::string>> surface_charge;
bool contains_id(const std::map<int, std::vector<std::string>> &map,
int id) const {
#if __cplusplus >= 202002L
return map.contains(id);
#else
return map.find(id) != map.end();
#endif
}
std::set<std::string> surface_primaries; std::set<std::string> surface_primaries;
using RawMap = std::map<int, std::pair<essential_names, std::vector<double>>>; using RawMap = std::map<int, std::pair<essential_names, std::vector<double>>>;

View File

@ -4,7 +4,89 @@
#include <span> #include <span>
#include <sstream> #include <sstream>
void PhreeqcEngine::init_wrappers(const POETInitCell &cell) { #include <IPhreeqc.hpp>
#include <Phreeqc.h>
#include "Wrapper/EquilibriumWrapper.hpp"
#include "Wrapper/ExchangeWrapper.hpp"
#include "Wrapper/KineticWrapper.hpp"
#include "Wrapper/SolutionWrapper.hpp"
#include "Wrapper/SurfaceWrapper.hpp"
class PhreeqcEngine::Impl : public IPhreeqc {
public:
Impl(const POETConfig &config);
void init_wrappers(const POETInitCell &cell);
void run(double time_step);
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;
};
PhreeqcEngine::PhreeqcEngine(const POETConfig &config)
: impl(std::make_unique<Impl>(config)) {}
PhreeqcEngine::Impl::Impl(const POETConfig &config) {
this->LoadDatabaseString(config.database.c_str());
this->RunString(config.input_script.c_str());
this->init_wrappers(config.cell);
}
PhreeqcEngine::~PhreeqcEngine() = default;
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->impl->set_essential_values(cell_data);
this->impl->run(time_step);
this->impl->get_essential_values(cell_data);
}
void PhreeqcEngine::Impl::run(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 1\n -time_step " + time_ss.str() + "\nEND\n";
this->RunString(runs_string.c_str());
}
void PhreeqcEngine::Impl::init_wrappers(const POETInitCell &cell) {
// Solutions // Solutions
this->solutionWrapperPtr = this->solutionWrapperPtr =
@ -41,7 +123,7 @@ void PhreeqcEngine::init_wrappers(const POETInitCell &cell) {
} }
} }
void PhreeqcEngine::get_essential_values(std::span<double> &data) { void PhreeqcEngine::Impl::get_essential_values(std::span<double> &data) {
this->solutionWrapperPtr->get(data); this->solutionWrapperPtr->get(data);
@ -76,54 +158,9 @@ void PhreeqcEngine::get_essential_values(std::span<double> &data) {
data.subspan(offset, this->surfaceWrapperPtr->size())}; data.subspan(offset, this->surfaceWrapperPtr->size())};
this->surfaceWrapperPtr->get(surf_span); 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) { void PhreeqcEngine::Impl::set_essential_values(const std::span<double> &data) {
this->solutionWrapperPtr->set(data); this->solutionWrapperPtr->set(data);
this->PhreeqcPtr->initial_solutions_poet(1); this->PhreeqcPtr->initial_solutions_poet(1);
@ -159,71 +196,4 @@ void PhreeqcEngine::set_essential_values(const std::span<double> &data) {
data.subspan(offset, this->surfaceWrapperPtr->size())}; data.subspan(offset, this->surfaceWrapperPtr->size())};
this->surfaceWrapperPtr->set(surf_span); 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(time_step);
this->get_essential_values(cell_data);
}
void PhreeqcEngine::run(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 1\n -time_step " + time_ss.str() + "\nEND\n";
this->RunString(runs_string.c_str());
} }

View File

@ -1,9 +1,9 @@
#include "EquilibriumWrapper.hpp"
#include "ExchangeWrapper.hpp"
#include "KineticWrapper.hpp"
#include "PhreeqcInit.hpp" #include "PhreeqcInit.hpp"
#include "SolutionWrapper.hpp" #include "Wrapper/EquilibriumWrapper.hpp"
#include "SurfaceWrapper.hpp" #include "Wrapper/ExchangeWrapper.hpp"
#include "Wrapper/KineticWrapper.hpp"
#include "Wrapper/SolutionWrapper.hpp"
#include "Wrapper/SurfaceWrapper.hpp"
#include "global_structures.h" #include "global_structures.h"
#include <cassert> #include <cassert>