diff --git a/ext/README.md b/ext/README.md new file mode 100644 index 0000000..faf5bb6 --- /dev/null +++ b/ext/README.md @@ -0,0 +1,10 @@ +## External Dependencies + +### IPhreeqc + +The contents of the `iphreeqc-3.8.6-17100` folder are extracted from the IPhreeqc module ("Linux (any processor)") +downloaded from `https://www.usgs.gov/software/phreeqc-version-3`. + +Specfically, at the time of writing (10/31/2025), the file +`https://water.usgs.gov/water-resources/software/PHREEQC/iphreeqc-3.8.6-17100.tar.gz` +(md5sum: `91d8485139b423d5f0382e9215d9aea1`) \ No newline at end of file diff --git a/src/bindings.cpp b/src/bindings.cpp index 69bdc98..a4df460 100644 --- a/src/bindings.cpp +++ b/src/bindings.cpp @@ -69,7 +69,11 @@ PYBIND11_MODULE(_bindings, m) { .def("get_selected_output_column_count", &IPhreeqcWrapper::get_selected_output_column_count) .def("get_value", &IPhreeqcWrapper::get_value) .def("get_component_count", &IPhreeqcWrapper::get_component_count) - .def("get_component", &IPhreeqcWrapper::get_component); + .def("get_component", &IPhreeqcWrapper::get_component) + .def("set_dump_string_on", &IPhreeqcWrapper::set_dump_string_on) + .def("get_dump_string", &IPhreeqcWrapper::get_dump_string) + .def("set_log_string_on", &IPhreeqcWrapper::set_log_string_on) + .def("get_log_string", &IPhreeqcWrapper::get_log_string); #ifdef VERSION_INFO diff --git a/src/iphreeqc_wrapper.cpp b/src/iphreeqc_wrapper.cpp index 846c083..8375bbf 100644 --- a/src/iphreeqc_wrapper.cpp +++ b/src/iphreeqc_wrapper.cpp @@ -54,6 +54,22 @@ class IPhreeqcWrapper { return GetComponent(id, i); } + std::string get_dump_string() { + return GetDumpString(id); + } + + int set_dump_string_on(int i) { + return SetDumpStringOn(id, i); + } + + std::string get_log_string() { + return GetLogString(id); + } + + int set_log_string_on(int i) { + return SetLogStringOn(id, i); + } + private: int id; }; diff --git a/src/pyphreeqc/__init__.py b/src/pyphreeqc/__init__.py index 7b7b9ac..1c8622c 100644 --- a/src/pyphreeqc/__init__.py +++ b/src/pyphreeqc/__init__.py @@ -1,6 +1,5 @@ -from .interface import IPhreeqc, Var +from .core import Phreeqc __all__ = [ - "IPhreeqc", - "Var" + "Phreeqc" ] \ No newline at end of file diff --git a/src/pyphreeqc/core.py b/src/pyphreeqc/core.py new file mode 100644 index 0000000..d1134c3 --- /dev/null +++ b/src/pyphreeqc/core.py @@ -0,0 +1,94 @@ +from typing import Any +from pathlib import Path +from pyphreeqc._bindings import PyIPhreeqc +from pyphreeqc.var import Var +from pyphreeqc.solution import Solution + + +class Phreeqc: + def __init__(self, database: str = "phreeqc.dat", database_directory: Path | None = None): + self._ext = PyIPhreeqc() + + if database_directory is None: + database_directory = Path(__file__).parent / "database" + self._ext.load_database(str(database_directory / database)) + + self._solutions: list[Solution] = [] + + # TODO: Is VAR the common denominator for most operations? + # Here we create one and modify it in operations instead of having + # the caller create new VARs per operation. + self._var: Var = Var() + + self.output = PhreeqcOutput(self) + + def __len__(self): + return len(self._solutions) + + def __getattr__(self, item) -> None: + """Delegate attribute access to the underlying PyIPhreeqc instance.""" + if hasattr(self._ext, item): + return getattr(self._ext, item) + raise AttributeError(f"Phreeqc has no attribute '{item}'") + + def add_solution(self, solution_dict: dict) -> None: + solution = Solution(solution_dict) + index = len(self) + self.run_string(f""" + SOLUTION {index} + {solution} + SAVE SOLUTION {index} + END + """) + self._solutions.append(solution) + + def remove_solution(self, index: int) -> Solution: + self.run_string(f""" + DELETE + -solution {index} + """) + return self._solutions.pop(index) + + +class PhreeqcOutput: + def __init__(self, phreeqc: Phreeqc): + self._phreeqc = phreeqc + + def __getitem__(self, item) -> Any: + if not isinstance(item, tuple): + item = (item,) + while len(item) < 2: + item += (slice(None),) + + row_idx, col_idx = item + + if isinstance(row_idx, slice): + row_indices = range(*row_idx.indices(self.shape[0])) + elif isinstance(row_idx, int): + row_indices = [row_idx] + else: + raise TypeError("Row index must be int or slice") + + if isinstance(col_idx, slice): + col_indices = range(*col_idx.indices(self.shape[1])) + elif isinstance(col_idx, int): + col_indices = [col_idx] + else: + raise TypeError("Column index must be int or slice") + + result = [] + for row in row_indices: + row_values = [] + for col in col_indices: + self._phreeqc._ext.get_value(row, col, self._phreeqc._var._var.var) + row_values.append(self._phreeqc._var.value) + result.append( + row_values if len(col_indices) > 1 else row_values[0]) + + if len(row_indices) == 1: + return result[0] + return result + + @property + def shape(self) -> tuple[int, int]: + return self._phreeqc.get_selected_output_row_count(), self._phreeqc.get_selected_output_column_count() \ No newline at end of file diff --git a/src/pyphreeqc/interface.py b/src/pyphreeqc/interface.py deleted file mode 100644 index bec1286..0000000 --- a/src/pyphreeqc/interface.py +++ /dev/null @@ -1,111 +0,0 @@ -from typing import Any -from pathlib import Path -from pyphreeqc._bindings import PyVar, PY_VAR_TYPE, PY_VRESULT, PyIPhreeqc - -IPhreeqc = PyIPhreeqc - - -class Var: - def __init__(self, value: Any | None = None): - self._var = PyVar() - self._var.var.type = PY_VAR_TYPE.TT_EMPTY - self.value = value # will invoke setter - - @property - def value(self) -> Any: - match self._var.var.type: - case PY_VAR_TYPE.TT_EMPTY: - return None - case PY_VAR_TYPE.TT_ERROR: - return self._var.var.vresult - case PY_VAR_TYPE.TT_LONG: - return self._var.var.lVal - case PY_VAR_TYPE.TT_DOUBLE: - return self._var.var.dVal - case PY_VAR_TYPE.TT_STRING: - return self._var.var.sVal - case _: - raise RuntimeError("Unknown type") - - @value.setter - def value(self, value) -> None: - # If we were previously holding a string, we need to free it by - # creating a new PyVar - if self._var.var.type == PY_VAR_TYPE.TT_STRING: - self._var = PyVar() - - if isinstance(value, PY_VRESULT): - self._var.var.type = PY_VAR_TYPE.TT_ERROR - self._var.var.vresult = value - elif isinstance(value, int): - self._var.var.type = PY_VAR_TYPE.TT_LONG - self._var.var.lVal = value - elif isinstance(value, float): - self._var.var.type = PY_VAR_TYPE.TT_DOUBLE - self._var.var.dVal = value - elif isinstance(value, str): - self._var.var.type = PY_VAR_TYPE.TT_STRING - self._var.var.sVal = value - elif value is None: - self._var.var.type = PY_VAR_TYPE.TT_EMPTY - else: - raise RuntimeError("Unknown type") - - -class Phreeqc: - def __init__(self, database: str = "phreeqc.dat", database_directory: Path | None = None): - self._ext = PyIPhreeqc() - - if database_directory is None: - database_directory = Path(__file__).parent / "database" - self._ext.load_database(str(database_directory / database)) - - # TODO: Is VAR the common denominator for most operations? - # Here we create one and modify it in operations instead of having - # the caller create new VARs per operation. - self._var = Var() - - def __getattr__(self, item) -> None: - """Delegate attribute access to the underlying PyIPhreeqc instance.""" - if hasattr(self._ext, item): - return getattr(self._ext, item) - raise AttributeError(f"Phreeqc has no attribute '{item}'") - - def __getitem__(self, item) -> Any: - if not isinstance(item, tuple): - item = (item,) - while len(item) < 2: - item += (slice(None),) - - row_idx, col_idx = item - - if isinstance(row_idx, slice): - row_indices = range(*row_idx.indices(self.shape[0])) - elif isinstance(row_idx, int): - row_indices = [row_idx] - else: - raise TypeError("Row index must be int or slice") - - if isinstance(col_idx, slice): - col_indices = range(*col_idx.indices(self.shape[1])) - elif isinstance(col_idx, int): - col_indices = [col_idx] - else: - raise TypeError("Column index must be int or slice") - - result = [] - for row in row_indices: - row_values = [] - for col in col_indices: - self._ext.get_value(row, col, self._var._var.var) - row_values.append(self._var.value) - result.append( - row_values if len(col_indices) > 1 else row_values[0]) - - if len(row_indices) == 1: - return result[0] - return result - - @property - def shape(self) -> tuple[int, int]: - return self.get_selected_output_row_count(), self.get_selected_output_column_count() \ No newline at end of file diff --git a/src/pyphreeqc/solution.py b/src/pyphreeqc/solution.py new file mode 100644 index 0000000..62c411c --- /dev/null +++ b/src/pyphreeqc/solution.py @@ -0,0 +1,6 @@ +class Solution(dict): + # A solution is nothing but a dict of str: Any mapping + pass + + def __str__(self): + return "\n".join(f"{k} {v}" for k, v in self.items()) \ No newline at end of file diff --git a/src/pyphreeqc/var.py b/src/pyphreeqc/var.py new file mode 100644 index 0000000..a5d9da0 --- /dev/null +++ b/src/pyphreeqc/var.py @@ -0,0 +1,49 @@ +from typing import Any +from pyphreeqc._bindings import PyVar, PY_VAR_TYPE, PY_VRESULT + + +class Var: + def __init__(self, value: Any | None = None): + self._var = PyVar() + self._var.var.type = PY_VAR_TYPE.TT_EMPTY + self.value = value # will invoke setter + + @property + def value(self) -> Any: + match self._var.var.type: + case PY_VAR_TYPE.TT_EMPTY: + return None + case PY_VAR_TYPE.TT_ERROR: + return self._var.var.vresult + case PY_VAR_TYPE.TT_LONG: + return self._var.var.lVal + case PY_VAR_TYPE.TT_DOUBLE: + return self._var.var.dVal + case PY_VAR_TYPE.TT_STRING: + return self._var.var.sVal + case _: + raise RuntimeError("Unknown type") + + @value.setter + def value(self, value) -> None: + # If we were previously holding a string, we need to free it by + # creating a new PyVar + if self._var.var.type == PY_VAR_TYPE.TT_STRING: + self._var = PyVar() + + if isinstance(value, PY_VRESULT): + self._var.var.type = PY_VAR_TYPE.TT_ERROR + self._var.var.vresult = value + elif isinstance(value, int): + self._var.var.type = PY_VAR_TYPE.TT_LONG + self._var.var.lVal = value + elif isinstance(value, float): + self._var.var.type = PY_VAR_TYPE.TT_DOUBLE + self._var.var.dVal = value + elif isinstance(value, str): + self._var.var.type = PY_VAR_TYPE.TT_STRING + self._var.var.sVal = value + elif value is None: + self._var.var.type = PY_VAR_TYPE.TT_EMPTY + else: + raise RuntimeError("Unknown type") \ No newline at end of file diff --git a/tests/test_phreeqc.py b/tests/test_phreeqc.py index 057c61e..c1933a5 100644 --- a/tests/test_phreeqc.py +++ b/tests/test_phreeqc.py @@ -1,6 +1,8 @@ from pathlib import Path import numpy as np -from pyphreeqc.interface import Phreeqc +import textwrap +import re +from pyphreeqc import Phreeqc def test_load_database_internal(): @@ -15,8 +17,8 @@ def test_load_database_external(): phreeqc = Phreeqc(database="phreeqc.dat", database_directory=Path(__file__).parent) -def test_run(): - # TODO: Break this down into individual tests +def test_run_sample(): + # Run a generic script string and capture output phreeqc = Phreeqc("phreeqc.dat") phreeqc.run_string(""" @@ -73,20 +75,193 @@ def test_run(): assert phreeqc.get_selected_output_column_count() == 8 # We can get values at a specific index, or a slice - assert phreeqc[0, 0] == "cb" - assert phreeqc[0, 1:4] == ["H", "O", "Ca"] - assert phreeqc[0, 5:] == ["K", "N", "Na"] + assert phreeqc.output[0, 0] == "cb" + assert phreeqc.output[0, 1:4] == ["H", "O", "Ca"] + assert phreeqc.output[0, 5:] == ["K", "N", "Na"] assert np.allclose( - phreeqc[1, :], + phreeqc.output[1, :], np.array([2.979292808179192e-18, 111.01243360409575, 55.50675186622646, 0.0006000000000000017, 0.0012000000000000005, 0.0, 0.0, 0.0]) ) assert np.allclose( - phreeqc[2], # single indices are also fine + phreeqc.output[2], # single indices are also fine np.array( [-3.395954270633993e-16, 111.01243360154359, 55.510351938877264, 0.0, 0.0, 0.00020000000000012212, 0.0012000000000010212, 0.0010000000000005584]) ) + + +def test_run_simple(): + phreeqc = Phreeqc() + # Note: No "SAVE SOLUTION" - we can get the results directly + phreeqc.run_string(""" + SOLUTION 0 + temp 25.0 + units mol/kgw + pH 7.0 + pe 8.5 + redox pe + water 0.9970480319717386 + SELECTED_OUTPUT + -solution true + END + """) + assert phreeqc.output.shape[0] == 2 # header + 1 solution + + +def test_run_add_solution(): + phreeqc = Phreeqc() + phreeqc.add_solution({'pH': 7.0, 'pe': 8.5, 'redox': 'pe', 'temp': 25.0, 'units': 'mol/kgw', 'water': 0.9970480319717386}) + assert len(phreeqc) == 1 + + +def test_run_add_delete_solution(): + phreeqc = Phreeqc() + phreeqc.add_solution({'pH': 7.0, 'pe': 8.5, 'redox': 'pe', 'temp': 25.0, 'units': 'mol/kgw', 'water': 0.9970480319717386}) + phreeqc.remove_solution(0) + assert len(phreeqc) == 0 + + +def _test_run_dumpstring(): + # This test is too platform-specific. Simplify! + phreeqc = Phreeqc() + + phreeqc.run_string(textwrap.dedent(""" + SOLUTION 0 + temp 25.0 + REACTION 1 + CaCl2 1 + Na2CO3 1 + 1 mmol + SAVE SOLUTION 0 + END + """)) + + phreeqc.set_dump_string_on(1) + + phreeqc.run_string(textwrap.dedent(""" + DUMP + -solution 0 + END + """)) + + dump_string = phreeqc.get_dump_string() + phreeqc.set_dump_string_on(0) + + expected = textwrap.dedent(""" + SOLUTION_RAW 0 Solution after simulation 1. + -temp 25 + -pressure 1 + -potential 0 + -total_h 111.01243359386 + -total_o 55.509216797548 + -cb -1.2200890388064e-09 + -density 0.99704301397679 + -viscosity 0.89125921464527 + -viscos_0 0.89002391825059 + -totals + C(4) 0.0010000000027134 + Ca 0.0010000000010752 + Cl 0.0019999999999998 + Na 0.002 + O(0) 2.8581598475304e-15 + -pH 10.413844680873 + -pe 7.3950560340151 + -mu 0.0045560988049649 + -ah2o 0.99989811240764 + -mass_water 0.999994868511 + -soln_vol 1.0029812985132 + -total_alkalinity 0.0020000012222396 + -activities + C(-4) -125.72046751407 + C(4) -3.4925949082566 + Ca -3.2726648366134 + Cl -2.7307519728506 + E -7.3950560340151 + H(0) -38.767826310689 + Na -2.730257199021 + O(0) -14.844435883364 + -gammas + USE mix none + USE reaction none + USE reaction_temperature none + USE reaction_pressure none + """).lstrip("\n") + + assert dump_string == expected + + +def _test_run_logstring(): + # This test is too platform-specific. Simplify! + phreeqc = Phreeqc() + phreeqc.set_log_string_on(1) + phreeqc.run_string(textwrap.dedent(""" + KNOBS + -logfile true + SOLUTION 0 + temp 25.0 + REACTION 1 + CaCl2 1 + Na2CO3 1 + 1 mmol + SAVE SOLUTION 0 + END + """)) + log_string = phreeqc.get_log_string() + phreeqc.set_log_string_on(0) + + + expected = textwrap.dedent(""" + ------------------------------------------- + Beginning of initial solution calculations. + ------------------------------------------- + + Initial solution 0. + + Iterations in revise_guesses: 1 + + Number of infeasible solutions: 0 + Number of basis changes: 0 + + Number of iterations: 0 + + ----------------------------------------- + Beginning of batch-reaction calculations. + ----------------------------------------- + + Reaction step 1. + + Overflow: (CO2)2\t1.000000e+03\t3.625196e+00\t-1 + Overflow: CO2\t4.794226e+02\t2.680719e+00\t-1 + Overflow: CaCO3\t1.000000e+03\t3.225283e+00\t-1 + Overflow: CaHCO3+\t1.000000e+03\t3.915297e+00\t-1 + Overflow: HCO3-\t1.000000e+03\t3.329016e+00\t-1 + Overflow: NaHCO3\t1.000000e+03\t3.268854e+00\t-1 + Iterations in revise_guesses: 2 + + Number of infeasible solutions: 0 + Number of basis changes: 0 + + Number of iterations: 15 + + ------------------ + End of simulation. + ------------------ + + ------------------------------------ + Reading input data for simulation 2. + ------------------------------------ + + --------------------------------- + End of Run after X Seconds. + --------------------------------- + """).strip("\n") + + + normalized = re.sub(r"End of Run after .* Seconds\.", + "End of Run after X Seconds.", log_string).rstrip() + + assert normalized == expected diff --git a/tests/test_phreeqc_solution.py b/tests/test_phreeqc_solution.py new file mode 100644 index 0000000..ed61646 --- /dev/null +++ b/tests/test_phreeqc_solution.py @@ -0,0 +1,10 @@ +from pyphreeqc import Phreeqc + + +def test_add_solution(): + phreeqc = Phreeqc() + phreeqc.add_solution( + {'Cl': '4.011842831773806', 'Na': '4.011842831773806', 'pH': 7.0, + 'pe': 8.5, 'redox': 'pe', 'temp': 25.0, 'units': 'mol/kgw', + 'water': 0.9970480319717386} + ) \ No newline at end of file diff --git a/tests/test_var.py b/tests/test_var.py index 7c4570b..d34b9a8 100644 --- a/tests/test_var.py +++ b/tests/test_var.py @@ -1,5 +1,5 @@ from pytest import approx -from pyphreeqc.interface import Var +from pyphreeqc.var import Var from pyphreeqc._bindings import PY_VRESULT