max/iphreeqc
1
0
Fork 0
mirror of https://git.gfz-potsdam.de/naaice/iphreeqc.git synced 2025-12-15 16:18:22 +01:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'ml/add_h0_o0' into 'poet'

Browse Source 
feat: add H(0) and O(0) to solution list

See merge request naaice/iphreeqc!11
This commit is contained in:
Max Lübke 2024-10-23 13:08:28 +02:00
commit 017d6bd69d
21 changed files with 1795 additions and 1818 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
.gitlab-ci.yml
View File

@ -2,7 +2,7 @@ image: ubuntu:24.04
before_script:
- apt-get update -y
- apt-get install -y build-essential cmake git
- apt-get install -y build-essential cmake git ninja-build
stages:
- test
@ -11,8 +11,8 @@ test:
stage: test
script:
- mkdir _build && cd _build
- cmake -DBUILD_SHARED_LIBS=ON -DCMAKE_BUILD_TYPE=Release ..
- make -j$(nproc)
- cmake -DBUILD_SHARED_LIBS=ON -DCMAKE_BUILD_TYPE=Release -G Ninja ..
- ninja
- ctest --output-junit test_results.xml
artifacts:
when: always

2
CMakeLists.txt
View File

@ -447,7 +447,7 @@ if (STANDALONE_BUILD)
gtest_hide_internal_symbols
)
FetchContent_GetProperties(googletest)
FetchContent_MakeAvailable(googletest)
if (NOT googletest_POPULATED)
FetchContent_Populate(googletest)
add_subdirectory(${googletest_SOURCE_DIR} ${googletest_BINARY_DIR})

19
poet/include/PhreeqcMatrix.hpp
View File

@ -9,6 +9,7 @@
#include <vector>
#include "IPhreeqc.hpp"
#include "global_structures.h"
/**
* @brief Class for storing information from Phreeqc
@ -166,19 +167,18 @@ public:
};
/**
* @brief Get all found solution names
* @brief Get all found solution names.
*
* @return std::vector<std::string> Vector containing all solution names.
*/
std::vector<std::string>
getSolutionNames(bool include_h_o_charge = false) const;
std::vector<std::string> getSolutionNames() const;
/**
* @brief Get solution total names of all found solutions (excluding H, O,
* Charge)
* Charge, H(0), O(0))
*
* @return std::vector<std::string> Names of all found solutions (excluding H,
* O, Charge)
* O, Charge, H(0), O(0))
*/
std::vector<std::string> getSolutionPrimaries() const;
@ -259,6 +259,15 @@ public:
*/
std::string getDatabase() const;
/**
* @brief Check if a cell with given ID exists in the PhreeqcMatrix.
*
* @param cell_id ID of the cell (user id from Phreeqc script) to check for.
* @return true Entry exists
* @return false Entry doesn't exist
*/
bool checkIfExists(int cell_id) const;
private:
std::map<int, std::vector<element>> _m_map;
std::map<int, std::vector<base_names>> _m_internal_names;

12
poet/src/Engine.cpp
View File

@ -80,6 +80,11 @@ replaceRawKeywordID(const std::string &raw_dump_string) {
}
PhreeqcEngine::Impl::Impl(const PhreeqcMatrix &pqc_mat, const int cell_id) {
if (!pqc_mat.checkIfExists(cell_id)) {
throw std::invalid_argument("Cell ID does not exist in PhreeqcMatrix");
}
this->LoadDatabaseString(pqc_mat.getDatabase().c_str());
const std::string pqc_string =
@ -100,6 +105,11 @@ PhreeqcEngine::Impl::Impl(const PhreeqcMatrix &pqc_mat, const int cell_id) {
void PhreeqcEngine::runCell(std::vector<double> &cell_values,
double time_step) {
if (time_step < 0) {
throw std::invalid_argument("Time step must be positive");
}
// skip ID
std::span<double> cell_data{cell_values.begin() + 1, cell_values.end()};
@ -194,7 +204,7 @@ void PhreeqcEngine::Impl::get_essential_values(std::span<double> &data) {
void PhreeqcEngine::Impl::set_essential_values(const std::span<double> &data) {
this->solutionWrapperPtr->set(data);
this->PhreeqcPtr->initial_solutions_poet(1);
// this->PhreeqcPtr->initial_solutions_poet(1);
std::size_t offset = this->solutionWrapperPtr->size();

15
poet/src/PhreeqcMatrix/Access.cpp
View File

@ -157,25 +157,18 @@ PhreeqcMatrix::STLExport PhreeqcMatrix::get(VectorExportType type,
return result;
}
std::vector<std::string>
PhreeqcMatrix::getSolutionNames(bool include_h_o_charge) const {
std::vector<std::string> PhreeqcMatrix::getSolutionNames() const {
std::vector<std::string> names;
if (include_h_o_charge) {
names.push_back("H");
names.push_back("O");
names.push_back("Charge");
}
const auto &first_element = _m_map.begin()->second;
for (std::size_t i = 3; i < _m_map.begin()->second.size(); i++) {
for (const auto &element : _m_map.begin()->second) {
// assuming the element vector always starts with the solution components
if (first_element[i].type != PhreeqcComponent::SOLUTION) {
if (element.type != PhreeqcComponent::SOLUTION) {
break;
}
names.push_back(first_element[i].name);
names.push_back(element.name);
}
return names;

7
poet/src/PhreeqcMatrix/Init.cpp
View File

@ -17,14 +17,12 @@
#include <utility>
#include <vector>
constexpr std::size_t SKIP_H_O_CB = 3;
static std::vector<std::string> dump_solution_names(cxxSolution *solution,
Phreeqc *phreeqc) {
std::vector<std::string> placeholder;
return phreeqc->find_all_valence_states(
SolutionWrapper::names(solution, placeholder), SKIP_H_O_CB);
SolutionWrapper::names(solution, placeholder));
}
template <enum PhreeqcMatrix::PhreeqcComponent comp, class T>
@ -122,8 +120,7 @@ create_vector_from_phreeqc(Phreeqc *phreeqc, int id,
// Solution
SolutionWrapper sol_wrapper(
Utilities::Rxn_find(phreeqc->Get_Rxn_solution_map(), id),
{solution_names.begin() + SKIP_H_O_CB, solution_names.end()});
Utilities::Rxn_find(phreeqc->Get_Rxn_solution_map(), id), solution_names);
base_add_to_element_vector<PhreeqcMatrix::PhreeqcComponent::SOLUTION>(
sol_wrapper, solution_names, elements);

7
poet/src/PhreeqcMatrix/Misc.cpp
View File

@ -1,5 +1,4 @@
#include "PhreeqcMatrix.hpp"
#include <cstddef>
std::vector<int> PhreeqcMatrix::getIds() const {
std::vector<int> ids;
@ -50,4 +49,8 @@ std::string PhreeqcMatrix::getDumpStringsPQI(int cell_id) const {
return dump_string;
}
std::string PhreeqcMatrix::getDatabase() const { return _m_database; }
std::string PhreeqcMatrix::getDatabase() const { return _m_database; }
bool PhreeqcMatrix::checkIfExists(int cell_id) const {
return this->_m_map.find(cell_id) != this->_m_map.end();
}

42
poet/src/Wrapper/EquilibriumWrapper.cpp
View File

@ -6,7 +6,7 @@
EquilibriumWrapper::EquilibriumWrapper(
cxxPPassemblage *ppassemblage,
const std::vector<std::string> &ppassemblage_names)
: ppassemblage(ppassemblage) {
: ppassemblage(ppassemblage), ppassemblage_order(ppassemblage_names) {
for (const auto &ppassemblage_name : ppassemblage_names) {
auto it = std::find_if(
ppassemblage->Get_pp_assemblage_comps().begin(),
@ -28,25 +28,49 @@ EquilibriumWrapper::EquilibriumWrapper(
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());
for (const auto &comp_name : ppassemblage_order) {
auto it = std::find_if(
equilibrium_comps.begin(), equilibrium_comps.end(),
[&](const auto &comp) { return comp->getCompName() == comp_name; });
comp->get(comp_span);
std::span<LDBLE> comp_span = data.subspan(offset, (*it)->size());
offset += comp->size();
(*it)->get(comp_span);
offset += (*it)->size();
}
// 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());
for (const auto &comp_name : ppassemblage_order) {
auto it = std::find_if(
equilibrium_comps.begin(), equilibrium_comps.end(),
[&](const auto &comp) { return comp->getCompName() == comp_name; });
comp->set(comp_span);
std::span<LDBLE> comp_span = data.subspan(offset, (*it)->size());
offset += comp->size();
(*it)->set(comp_span);
offset += (*it)->size();
}
// 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>

4
poet/src/Wrapper/EquilibriumWrapper.hpp
View File

@ -3,6 +3,7 @@
#include "PPassemblage.h"
#include "WrapperBase.hpp"
#include <memory>
#include <vector>
class EquilibriumWrapper : public WrapperBase {
public:
@ -30,9 +31,12 @@ private:
static std::vector<std::string> names(const cxxPPassemblageComp &comp);
std::string getCompName() const { return this->comp.Get_name(); }
private:
cxxPPassemblageComp &comp;
};
std::vector<std::unique_ptr<EquilibriumCompWrapper>> equilibrium_comps;
const std::vector<std::string> ppassemblage_order;
};

42
poet/src/Wrapper/SolutionWrapper.cpp
View File

@ -1,11 +1,13 @@
#include "SolutionWrapper.hpp"
#include "NameDouble.h"
#include <set>
#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;
cxxSolution *soln, const std::vector<std::string> &_solution_order)
: solution(soln), solution_order(_solution_order.begin() + NUM_ESSENTIALS,
_solution_order.end()) {
this->num_elements = _solution_order.size();
auto &totals = solution->Get_totals();
}
@ -14,6 +16,8 @@ void SolutionWrapper::get(std::span<LDBLE> &data) const {
data[0] = solution->Get_total_h();
data[1] = solution->Get_total_o();
data[2] = solution->Get_cb();
data[3] = solution->Get_total("H(0)");
data[4] = solution->Get_total("O(0)");
std::size_t i = NUM_ESSENTIALS;
for (const auto &tot_name : solution_order) {
@ -22,13 +26,21 @@ void SolutionWrapper::get(std::span<LDBLE> &data) const {
data[i++] = 0.0;
continue;
}
data[i++] = it->second;
data[i++] = it->second > 1e-25 ? it->second : 0.;
}
}
void SolutionWrapper::set(const std::span<LDBLE> &data) {
std::size_t i = NUM_ESSENTIALS;
cxxNameDouble new_totals;
const double &total_h = data[0];
const double &total_o = data[1];
const double &cb = data[2];
new_totals["H(0)"] = data[3];
new_totals["O(0)"] = data[4];
for (const auto &tot_name : solution_order) {
const double value = data[i++];
@ -38,7 +50,7 @@ void SolutionWrapper::set(const std::span<LDBLE> &data) {
new_totals[tot_name] = value;
}
this->solution->Update(data[0], data[1], data[2], new_totals);
this->solution->Update(total_h, total_o, cb, new_totals);
}
std::vector<std::string>
@ -46,18 +58,18 @@ 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");
names.insert(names.end(), ESSENTIALS.begin(), ESSENTIALS.end());
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);
std::set<std::string> names_set;
for (const auto &name : solution->Get_totals()) {
names_set.insert(name.first);
}
names.insert(names.end(), solution_order.begin(), solution_order.end());
for (const auto &to_erase : ESSENTIALS) {
// don't care if the element was not found
names_set.erase(to_erase);
}
names.insert(names.end(), names_set.begin(), names_set.end());
return names;
}

9
poet/src/Wrapper/SolutionWrapper.hpp
View File

@ -3,7 +3,9 @@
#include "NameDouble.h"
#include "Solution.h"
#include "WrapperBase.hpp"
#include <array>
#include <cstddef>
#include <string>
#include <vector>
class SolutionWrapper : public WrapperBase {
@ -18,9 +20,14 @@ public:
static std::vector<std::string>
names(cxxSolution *solution, std::vector<std::string> &solution_order);
std::vector<std::string> getEssentials() const;
private:
cxxSolution *solution;
const std::vector<std::string> solution_order;
static constexpr std::size_t NUM_ESSENTIALS = 3;
static constexpr std::array<std::string, 5> ESSENTIALS = {"H", "O", "Charge",
"H(0)", "O(0)"};
static constexpr std::size_t NUM_ESSENTIALS = ESSENTIALS.size();
};

10
poet/test/CMakeLists.txt
View File

@ -2,7 +2,9 @@ enable_testing()
add_executable(
poet_test
testPhreeqcEngine.cpp
testPhreeqcMatrix.cpp
utils.cpp
)
target_link_libraries(
@ -14,11 +16,11 @@ target_link_libraries(
target_include_directories(poet_test PRIVATE ${CMAKE_CURRENT_BINARY_DIR})
# read file and store in variable
file(READ "${PROJECT_SOURCE_DIR}/database/phreeqc.dat" POET_PHREEQCDAT_DB)
file(READ "${CMAKE_CURRENT_SOURCE_DIR}/barite_db.dat" POET_BARITE_DB)
file(READ "${CMAKE_CURRENT_SOURCE_DIR}/barite_het.pqi" POET_BARITE_PQI)
file(REAL_PATH "${PROJECT_SOURCE_DIR}/database/phreeqc.dat" POET_PHREEQCDAT_DB)
file(REAL_PATH "${CMAKE_CURRENT_SOURCE_DIR}/barite_db.dat" POET_BARITE_DB)
file(REAL_PATH "${CMAKE_CURRENT_SOURCE_DIR}/barite_het.pqi" POET_BARITE_PQI)
configure_file("${CMAKE_CURRENT_SOURCE_DIR}/testPhreeqcMatrix.hpp.in" "${CMAKE_CURRENT_BINARY_DIR}/testPhreeqcMatrix.hpp")
configure_file("${CMAKE_CURRENT_SOURCE_DIR}/testInput.hpp.in" "${CMAKE_CURRENT_BINARY_DIR}/testInput.hpp")
include(GoogleTest)
gtest_discover_tests(poet_test)

128
poet/test/testInput.hpp.in Normal file
View File

@ -0,0 +1,128 @@
#pragma once
#include <limits>
#include <string>
#include <vector>
#define POET_TEST(name) TEST(TestPOET, name)
namespace base_test {
const std::string script = 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)";
const std::vector<std::string> expected_names = {
"ID", "H", "O", "Charge", "H(0)", "O(0)",
"C(-4)", "C(4)", "Ca", "Cl", "Mg", "Na",
"Calcite_eq", "Calcite_si", "Dolomite_eq", "Dolomite_si"};
const std::vector<double> expected_values = {1,
111.01243522078478,
55.506323405120348,
-4.7919630342304069e-13,
0,
2.0180693312370871e-15,
0,
3.500625006800175e-05,
0.12131244646561848,
0.49999844804496646,
0.078722559784449683,
0.099999999999999978,
0.95741011331883141,
0,
0.021277440215550378,
0};
const std::vector<double> expected_errors = {
0, 1e-3, 1e-3, 1e-15, 1e-15, 1e-15, 1e-5, 1e-5,
1e-5, 1e-5, 1e-5, 1e-5, 1e-5, 0, 1e-5, 0};
const std::string phreeqc_database = R"database(@POET_PHREEQCDAT_DB@)database";
} // namespace base_test
namespace barite_test {
const std::string script = R"barite(@POET_BARITE_PQI@)barite";
const std::string database = R"barite(@POET_BARITE_DB@)barite";
const std::vector<std::string> expected_names = {
"ID", "H", "O", "Charge", "H(0)",
"O(0)", "Ba", "Cl", "S(-2)", "S(6)",
"Sr", "Barite", "Barite_p1", "Celestite", "Celestite_p1",
"Celestite_eq", "Celestite_si"};
const std::vector<double> expected_values_line_one = {
1,
111.01243359383071,
55.508698688362124,
-1.2153078399577636e-09,
0,
8.6371063688066983e-15,
1.0000001848805677e-12,
1.0000000116187218e-12,
0,
0.00062047270964839664,
0.00062047270964840271,
std::numeric_limits<double>::quiet_NaN(),
std::numeric_limits<double>::quiet_NaN(),
std::numeric_limits<double>::quiet_NaN(),
std::numeric_limits<double>::quiet_NaN(),
0.99937952729135193,
0};
const std::vector<double> expected_errors = {
0,
1e-3,
1e-3,
1e-15,
1e-15,
1e-15,
1e-5,
1e-5,
1e-5,
1e-5,
1e-5,
std::numeric_limits<double>::quiet_NaN(),
std::numeric_limits<double>::quiet_NaN(),
std::numeric_limits<double>::quiet_NaN(),
std::numeric_limits<double>::quiet_NaN(),
1e-5,
0};
const std::vector<std::string> expected_names_erased = {
"ID", "H", "O", "Charge", "H(0)",
"O(0)", "Ba", "Cl", "S(-2)", "S(6)",
"Sr", "Barite", "Barite_p1", "Celestite", "Celestite_p1"};
const std::vector<std::string> expected_names_subset = {
"ID", "H", "O", "Charge", "H(0)", "O(0)", "Ba",
"Cl", "S(-2)", "S(6)", "Sr", "Celestite_eq", "Celestite_si"};
} // namespace barite_test
namespace test_engine {
const std::string script = 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)";
const std::string phreeqc_database = R"database(@POET_PHREEQCDAT_DB@)database";
} // namespace test_engine

41
poet/test/testPhreeqcEngine.cpp Normal file
View File

@ -0,0 +1,41 @@
#include <stdexcept>
#include <testInput.hpp>
#include <gtest/gtest.h>
#include "PhreeqcEngine.hpp"
#include "PhreeqcMatrix.hpp"
#include "utils.hpp"
const std::string test_database = readFile(base_test::phreeqc_database);
POET_TEST(PhreeqcEngineConstructor) {
PhreeqcMatrix pqc_mat(test_database, base_test::script);
EXPECT_NO_THROW(PhreeqcEngine(pqc_mat, 1));
EXPECT_THROW(PhreeqcEngine(pqc_mat, 2), std::invalid_argument);
}
POET_TEST(PhreeqcEngineStep) {
PhreeqcMatrix pqc_mat(test_database, base_test::script);
PhreeqcEngine engine(pqc_mat, 1);
std::vector<double> cell_values = pqc_mat.get().values;
EXPECT_NO_THROW(engine.runCell(cell_values, 0));
EXPECT_NO_THROW(engine.runCell(cell_values, 100));
for (std::size_t i = 0; i < cell_values.size(); ++i) {
// ignore H(0) and O(0)
if (i == 4 || i == 5) {
continue;
}
EXPECT_NEAR(cell_values[i], base_test::expected_values[i],
base_test::expected_errors[i]);
}
EXPECT_THROW(engine.runCell(cell_values, -1), std::invalid_argument);
}

54
poet/test/testPhreeqcMatrix.cpp
View File

@ -2,26 +2,29 @@
#include <gtest/gtest.h>
#include <vector>
#include <testPhreeqcMatrix.hpp>
#include <testInput.hpp>
#include "PhreeqcMatrix.hpp"
#include "utils.hpp"
#define POET_TEST(name) TEST(TestPOET, name)
const std::string base_db = readFile(base_test::phreeqc_database);
POET_TEST(PhreeqcInit) {
EXPECT_NO_THROW(
PhreeqcMatrix(base_test::phreeqc_database, base_test::script));
EXPECT_NO_THROW(PhreeqcMatrix(base_db, base_test::script));
}
POET_TEST(PhreeqcMatrixOneSolution) {
PhreeqcMatrix pqc_mat(base_test::phreeqc_database, base_test::script);
PhreeqcMatrix pqc_mat(base_db, base_test::script);
const auto ids = pqc_mat.getIds();
EXPECT_EQ(ids.size(), 1);
EXPECT_EQ(ids[0], 1);
EXPECT_TRUE(pqc_mat.checkIfExists(1));
EXPECT_FALSE(pqc_mat.checkIfExists(2));
PhreeqcMatrix::STLExport exported_init = pqc_mat.get();
// ID + H,O,Charge + 4 Solutions + 4 Equil incl. params
EXPECT_EQ(exported_init.names.size(), 12);
// ID + H,O,Charge,H(0),O(0) + 6 Solutions + 4 Equil incl. params
EXPECT_EQ(exported_init.names.size(), 16);
EXPECT_EQ(exported_init.names, base_test::expected_names);
for (std::size_t i = 0; i < exported_init.values.size(); ++i) {
@ -44,8 +47,20 @@ POET_TEST(PhreeqcMatrixOneSolution) {
EXPECT_EQ(equilibrium, expected_equilibrium);
}
POET_TEST(PhreeqcMatrixBracketOperator) {
PhreeqcMatrix pqc_mat(base_db, base_test::script);
EXPECT_NO_THROW(pqc_mat(1, "H"));
EXPECT_NEAR(pqc_mat(1, "H"), base_test::expected_values[1], 1e-5);
EXPECT_ANY_THROW(pqc_mat(1, "J"));
EXPECT_ANY_THROW(pqc_mat(2, "H"));
}
const std::string barite_db = readFile(barite_test::database);
const std::string barite_script = readFile(barite_test::script);
POET_TEST(PhreeqcMatrixMultiSolution) {
PhreeqcMatrix pqc_mat(barite_test::database, barite_test::script);
PhreeqcMatrix pqc_mat(barite_db, barite_script);
const auto ids = pqc_mat.getIds();
EXPECT_EQ(ids.size(), 4);
@ -58,7 +73,7 @@ POET_TEST(PhreeqcMatrixMultiSolution) {
EXPECT_EQ(exported.names, barite_test::expected_names);
for (std::size_t i = 0; i < exported.names.size(); i++) {
if (i > 8 && i < 13) {
if (i > 10 && i < 15) {
EXPECT_TRUE(std::isnan(exported.values[i]));
continue;
}
@ -94,7 +109,7 @@ POET_TEST(PhreeqcMatrixMultiSolution) {
}
POET_TEST(PhreeqcMatrixCtor) {
PhreeqcMatrix pqc_mat(barite_test::database, barite_test::script);
PhreeqcMatrix pqc_mat(barite_db, barite_script);
PhreeqcMatrix pqc_mat_copy(pqc_mat);
PhreeqcMatrix pqc_mat_move(std::move(pqc_mat_copy));
@ -109,7 +124,7 @@ POET_TEST(PhreeqcMatrixCtor) {
EXPECT_EQ(exported.names, barite_test::expected_names);
for (std::size_t i = 0; i < exported.names.size(); i++) {
if (i > 8 && i < 13) {
if (i > 10 && i < 15) {
EXPECT_TRUE(std::isnan(exported.values[i]));
continue;
}
@ -119,7 +134,7 @@ POET_TEST(PhreeqcMatrixCtor) {
}
POET_TEST(PhreeqcMatrixOperator) {
PhreeqcMatrix pqc_mat(barite_test::database, barite_test::script);
PhreeqcMatrix pqc_mat(barite_db, barite_script);
PhreeqcMatrix pqc_mat_copy = pqc_mat;
const auto ids = pqc_mat_copy.getIds();
@ -134,7 +149,7 @@ POET_TEST(PhreeqcMatrixOperator) {
EXPECT_EQ(exported.names, barite_test::expected_names);
for (std::size_t i = 0; i < exported.names.size(); i++) {
if (i > 8 && i < 13) {
if (i > 10 && i < 15) {
EXPECT_TRUE(std::isnan(exported.values[i]));
continue;
}
@ -144,7 +159,7 @@ POET_TEST(PhreeqcMatrixOperator) {
}
POET_TEST(PhreeqcMatrixRvalueManipulation) {
PhreeqcMatrix pqc_mat(barite_test::database, barite_test::script);
PhreeqcMatrix pqc_mat(barite_db, barite_script);
PhreeqcMatrix pqc_erased = pqc_mat.erase({1});
@ -175,7 +190,7 @@ POET_TEST(PhreeqcMatrixRvalueManipulation) {
}
POET_TEST(PhreeqcMatrixColumnMajorExport) {
PhreeqcMatrix pqc_mat(barite_test::database, barite_test::script);
PhreeqcMatrix pqc_mat(barite_db, barite_script);
pqc_mat = pqc_mat.subset({2, 3});
@ -193,12 +208,3 @@ POET_TEST(PhreeqcMatrixColumnMajorExport) {
EXPECT_EQ(exported.values[0], 2);
EXPECT_EQ(exported.values[1], 3);
}
POET_TEST(PhreeqcMatrixBracketOperator) {
PhreeqcMatrix pqc_mat(base_test::phreeqc_database, base_test::script);
EXPECT_NO_THROW(pqc_mat(1, "H"));
EXPECT_NEAR(pqc_mat(1, "H"), base_test::expected_values[1], 1e-5);
EXPECT_ANY_THROW(pqc_mat(1, "J"));
EXPECT_ANY_THROW(pqc_mat(2, "H"));
}

105
poet/test/testPhreeqcMatrix.hpp.in
View File

@ -1,105 +0,0 @@
#pragma once
#include <limits>
#include <string>
#include <vector>
namespace base_test {
const std::string script = 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)";
const std::vector<std::string> expected_names = {
"ID", "H", "O", "Charge", "Ca", "Cl",
"Mg", "Na", "Calcite_eq", "Calcite_si", "Dolomite_eq", "Dolomite_si"};
const std::vector<double> expected_values = {1,
111.01243521533338,
55.506218386370165,
-4.7941959748097226e-13,
0.1,
0.5,
0.1,
0.1,
1,
0,
0,
0};
const std::vector<double> expected_errors = {
0, 1e-3, 1e-3, 1e-15, 1e-5, 1e-5, 1e-5, 1e-5, 1e-5, 0, 1e-5, 0};
const std::string phreeqc_database = R"database(@POET_PHREEQCDAT_DB@)database";
} // namespace base_test
namespace barite_test {
const std::string script = R"barite(@POET_BARITE_PQI@)barite";
const std::string database = R"barite(@POET_BARITE_DB@)barite";
const std::vector<std::string> expected_names = {"ID",
"H",
"O",
"Charge",
"Ba",
"Cl",
"S(-2)",
"S(6)",
"Sr",
"Barite",
"Barite_p1",
"Celestite",
"Celestite_p1",
"Celestite_eq",
"Celestite_si"};
const std::vector<double> expected_values_line_one = {
1,
111.01243359383071,
55.508698688362124,
-1.2153078399577636e-09,
1.0000001848805677e-12,
1.0000000116187218e-12,
0,
0.00062047270964839664,
0.00062047270964840271,
std::numeric_limits<double>::quiet_NaN(),
std::numeric_limits<double>::quiet_NaN(),
std::numeric_limits<double>::quiet_NaN(),
std::numeric_limits<double>::quiet_NaN(),
0.99937952729135193,
0};
const std::vector<double> expected_errors = {
0,
1e-3,
1e-3,
1e-15,
1e-5,
1e-5,
1e-5,
1e-5,
1e-5,
std::numeric_limits<double>::quiet_NaN(),
std::numeric_limits<double>::quiet_NaN(),
std::numeric_limits<double>::quiet_NaN(),
std::numeric_limits<double>::quiet_NaN(),
1e-5,
0};
const std::vector<std::string> expected_names_erased = {
"ID", "H", "O", "Charge", "Ba", "Cl", "S(-2)",
"S(6)", "Sr", "Barite", "Barite_p1", "Celestite", "Celestite_p1"};
const std::vector<std::string> expected_names_subset = {
"ID", "H", "O", "Charge", "Ba", "Cl",
"S(-2)", "S(6)", "Sr", "Celestite_eq", "Celestite_si"};
} // namespace barite_test

25
poet/test/utils.cpp Normal file
View File

@ -0,0 +1,25 @@
#include "utils.hpp"
#include <fstream>
#include <linux/limits.h>
#include <sstream>
#include <string>
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();
}

5
poet/test/utils.hpp Normal file
View File

@ -0,0 +1,5 @@
#pragma once
#include <string>
std::string readFile(const std::string &path);

21
src/phreeqcpp/GFZ.cpp
View File

@ -1,14 +1,27 @@
#include "Phreeqc.h"
#include <set>
const std::set<std::string> to_ignore = {"H", "O", "Charge", "H(0)", "O(0)"};
std::vector<std::string> Phreeqc::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);
const std::vector<std::string> &&solution_names) {
std::vector<std::string> solution_with_valences;
solution_with_valences.reserve(solution_names.size());
// 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++) {
auto is_ignored = [](const std::string &in) {
return (to_ignore.find(in) != to_ignore.end());
};
for (std::size_t i = 0; i < solution_names.size(); i++) {
if (is_ignored(solution_names[i])) {
solution_with_valences.push_back(solution_names[i]);
continue;
}
const auto master_primary =
master_bsearch_primary(solution_names[i].c_str());

3
src/phreeqcpp/Phreeqc.h
View File

@ -1842,8 +1842,7 @@ protected:
public:
std::vector<std::string>
find_all_valence_states(const std::vector<std::string> &&solution_names,
const std::size_t offset);
find_all_valence_states(const std::vector<std::string> &&solution_names);
static const class const_iso iso_defaults[];
static const int count_iso_defaults;

3056
src/phreeqcpp/Solution.cxx
View File

File diff suppressed because it is too large Load Diff