From cbde10904ca944f8653bfd9021a87660236c4ba2 Mon Sep 17 00:00:00 2001 From: Darth Vader Date: Fri, 3 Jan 2025 19:08:20 +0000 Subject: [PATCH] Squashed 'src/phast/PhreeqcRM/' changes from 166e4032..a731364d a731364d 188 remove pybind (#189) git-subtree-dir: src/phast/PhreeqcRM git-subtree-split: a731364d2ed32eb29079e65ba247970db3548f60 --- .github/workflows/ci.yml | 8 +- .gitmodules | 4 - README.md | 5 + pybind/CMakeLists.txt | 230 ---- pybind/docstrings.h | 124 -- pybind/pybind.cpp | 2211 -------------------------------- pybind/pybind.cxx | 986 -------------- pybind/pybind11 | 1 - pybind/tests/AdvectBMI_py.yaml | 1163 ----------------- pybind/tests/advect.pqi | 31 - pybind/tests/constants.py | 19 - pybind/tests/phreeqc.dat | 1930 ---------------------------- pybind/tests/test_basic.py | 396 ------ pybind/tests/test_get_value.py | 117 -- pybind/tests/test_grid_info.py | 90 -- pybind/tests/test_irf.py | 95 -- pybind/tests/test_set_value.py | 51 - 17 files changed, 11 insertions(+), 7450 deletions(-) delete mode 100644 pybind/CMakeLists.txt delete mode 100644 pybind/docstrings.h delete mode 100644 pybind/pybind.cpp delete mode 100644 pybind/pybind.cxx delete mode 160000 pybind/pybind11 delete mode 100644 pybind/tests/AdvectBMI_py.yaml delete mode 100644 pybind/tests/advect.pqi delete mode 100644 pybind/tests/constants.py delete mode 100644 pybind/tests/phreeqc.dat delete mode 100644 pybind/tests/test_basic.py delete mode 100644 pybind/tests/test_get_value.py delete mode 100644 pybind/tests/test_grid_info.py delete mode 100644 pybind/tests/test_irf.py delete mode 100644 pybind/tests/test_set_value.py diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 68cc1be13..4d7777848 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -130,7 +130,11 @@ jobs: valgrind: # if: ${{ false }} name: valgrind - runs-on: ubuntu-latest + strategy: + fail-fast: false + matrix: + os: [ubuntu-24.04, ubuntu-22.04] + runs-on: ${{ matrix.os }} steps: - uses: actions/checkout@v4 @@ -349,7 +353,7 @@ jobs: strategy: fail-fast: false matrix: - os: [macos-13, ubuntu-20.04, ubuntu-22.04, windows-2019, windows-2022] + os: [macos-13, ubuntu-20.04, ubuntu-22.04, ubuntu-24.04, windows-2019, windows-2022] build_shared_libs: [OFF, ON] phreeqcrm_with_yaml_cpp: [OFF, ON] BUILD_TYPE: [Debug, Release] diff --git a/.gitmodules b/.gitmodules index f2cd90866..e69de29bb 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,4 +0,0 @@ -[submodule "pybind/pybind11"] - path = pybind/pybind11 - url = https://github.com/pybind/pybind11.git - branch = stable diff --git a/README.md b/README.md index 79386211e..0d2459031 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,10 @@ # PhreeqcRM: A reaction module for transport simulators based on the geochemical model PHREEQC +![PyPI - Version](https://img.shields.io/pypi/v/phreeqcrm) +![Conda Version](https://img.shields.io/conda/v/conda-forge/phreeqcrm) +![PyPI - Python Version](https://img.shields.io/pypi/pyversions/phreeqcrm) + + Reaction module for reactive-transport simulators. Based on IPhreeqc, allows access to all PHREEQC reaction capabilities. Contains methods for initial and boundary conditions, running reactions on all model cells, transfer of data to and from the module, and parallelization by MPI or OpenMP. This is the development repository for PhreeqcRM. The official USGS distribution is available at [USGS Release Page](https://www.usgs.gov/software/phreeqc-version-3). diff --git a/pybind/CMakeLists.txt b/pybind/CMakeLists.txt deleted file mode 100644 index c78d12213..000000000 --- a/pybind/CMakeLists.txt +++ /dev/null @@ -1,230 +0,0 @@ -cmake_minimum_required(VERSION 3.4...3.18) - -project(pybindtest LANGUAGES CXX C) - -##set_property(GLOBAL PROPERTY PROJECT_LABEL "pybind11") - -set(PYTHON_TARGET_NAME phreeqcrm) - -set(PhreeqcRM_SOURCES - ../src/BMI_interface_F.cpp - ../src/BMI_interface_F.h - ../src/BMIVariant.cpp - ../src/BMIVariant.h - ../src/bmi.hxx - ../src/BMIPhreeqcRM.cpp - ../src/BMIPhreeqcRM.h - ../src/IPhreeqcPhast/IPhreeqc/CSelectedOutput.cpp - ../src/IPhreeqcPhast/IPhreeqc/CSelectedOutput.hxx - ../src/IPhreeqcPhast/IPhreeqc/CVar.hxx - ../src/IPhreeqcPhast/IPhreeqc/Debug.h - ../src/IPhreeqcPhast/IPhreeqc/ErrorReporter.hxx - ../src/IPhreeqcPhast/IPhreeqc/IPhreeqc_interface_F.cpp - ../src/IPhreeqcPhast/IPhreeqc/IPhreeqc_interface_F.h - ../src/IPhreeqcPhast/IPhreeqc/IPhreeqc.cpp - ../src/IPhreeqcPhast/IPhreeqc/IPhreeqc.h - ../src/IPhreeqcPhast/IPhreeqc/IPhreeqc.hpp - ../src/IPhreeqcPhast/IPhreeqc/IPhreeqcCallbacks.h - ../src/IPhreeqcPhast/IPhreeqc/IPhreeqcLib.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/advection.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/basicsubs.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/cl1.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/common/Parser.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/common/Parser.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/common/PHRQ_base.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/common/PHRQ_base.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/common/PHRQ_exports.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/common/PHRQ_io.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/common/PHRQ_io.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/common/phrqtype.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/common/Utils.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/common/Utils.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/cvdense.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/cvdense.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/cvode.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/cvode.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/cxxKinetics.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/cxxKinetics.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/cxxMix.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/cxxMix.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/dense.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/dense.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/Dictionary.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/Dictionary.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/dumper.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/dumper.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/Exchange.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/Exchange.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/ExchComp.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/ExchComp.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/GasComp.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/GasComp.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/gases.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/GasPhase.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/GasPhase.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/global_structures.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/input.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/integrate.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/inverse.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/ISolution.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/ISolution.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/ISolutionComp.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/ISolutionComp.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/isotopes.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/kinetics.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/KineticsComp.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/KineticsComp.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/mainsubs.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/model.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/NA.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/NameDouble.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/NameDouble.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/NumKeyword.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/NumKeyword.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/nvector_serial.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/nvector_serial.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/nvector.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/nvector.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/parse.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/PBasic.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/PBasic.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/phqalloc.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/phqalloc.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/Phreeqc.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/Phreeqc.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/PhreeqcKeywords/Keywords.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/PhreeqcKeywords/Keywords.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/PHRQ_io_output.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/pitzer_structures.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/pitzer.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/PPassemblage.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/PPassemblage.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/PPassemblageComp.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/PPassemblageComp.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/prep.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/Pressure.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/Pressure.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/print.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/Reaction.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/Reaction.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/read.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/ReadClass.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/readtr.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/runner.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/runner.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/SelectedOutput.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/SelectedOutput.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/Serializer.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/Serializer.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/sit.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/smalldense.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/smalldense.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/Solution.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/Solution.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/SolutionIsotope.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/SolutionIsotope.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/spread.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/SS.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/SS.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/SSassemblage.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/SSassemblage.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/SScomp.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/SScomp.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/step.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/StorageBin.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/StorageBin.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/StorageBinList.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/StorageBinList.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/structures.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/sundialsmath.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/sundialsmath.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/sundialstypes.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/Surface.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/Surface.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/SurfaceCharge.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/SurfaceCharge.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/SurfaceComp.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/SurfaceComp.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/System.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/System.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/tally.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/Temperature.cxx - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/Temperature.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/tidy.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/transport.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/Use.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/Use.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/UserPunch.cpp - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/UserPunch.h - ../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/utilities.cpp - ../src/IPhreeqcPhast/IPhreeqc/thread.h - ../src/IPhreeqcPhast/IPhreeqc/Var.c - ../src/IPhreeqcPhast/IPhreeqc/Var.h - ../src/IPhreeqcPhast/IPhreeqc/Version.h - ../src/IPhreeqcPhast/IPhreeqcPhast.cxx - ../src/IPhreeqcPhast/IPhreeqcPhast.h - ../src/IPhreeqcPhast/IPhreeqcPhastLib.cpp - ../src/IPhreeqcPhast/IPhreeqcPhastLib.h - ../src/IrmResult.h - ../src/PhreeqcRM.cpp - ../src/PhreeqcRM.h - ../src/BMI_interface_F.cpp - ../src/RM_interface_C.cpp - ../src/RM_interface_C.h - ../src/RM_interface_F.cpp - ../src/RM_interface_F.h - ../src/RMVARS.h - ../src/VarManager.cpp - ../src/VarManager.h - ../src/YAML_interface_F.cpp - ../src/YAML_interface_F.h - ../src/YAMLPhreeqcRM.cpp - ../src/YAMLPhreeqcRM.h -) - -# option(PHREEQCRM_WITH_YAML_CPP "Build with yaml-cpp support" OFF) -# if(PHREEQCRM_WITH_YAML_CPP) - find_package(yaml-cpp REQUIRED) -# endif() - -add_subdirectory(pybind11) - -pybind11_add_module(${PYTHON_TARGET_NAME} pybind.cxx docstrings.h ${PhreeqcRM_SOURCES}) -##pybind11_add_module(${PYTHON_TARGET_NAME} pybind.cpp) -##pybind11_add_module(cmake_example pybind.cpp) - -##set_target_properties(${PYTHON_TARGET_NAME} PROPERTIES TARGET_NAME "MyCustomTarget") - -# pybind11 defs -target_compile_definitions(${PYTHON_TARGET_NAME} PRIVATE WITH_PYBIND11) - -# iphreeqc defs -target_compile_definitions(${PYTHON_TARGET_NAME} PRIVATE SWIG_SHARED_OBJ) -target_compile_definitions(${PYTHON_TARGET_NAME} PRIVATE USE_PHRQ_ALLOC) -target_compile_definitions(${PYTHON_TARGET_NAME} PRIVATE IPhreeqc_EXPORTS) - - -# if(PHREEQCRM_WITH_YAML_CPP AND yaml-cpp_FOUND) - target_compile_definitions(${PYTHON_TARGET_NAME} PRIVATE USE_YAML) - target_link_libraries(${PYTHON_TARGET_NAME} PRIVATE yaml-cpp) -# endif() - -target_include_directories(${PYTHON_TARGET_NAME} - PRIVATE - "../src" - "../src/IPhreeqcPhast" - "../src/IPhreeqcPhast/IPhreeqc" - "../src/IPhreeqcPhast/IPhreeqc/phreeqcpp" - "../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/common" - "../src/IPhreeqcPhast/IPhreeqc/phreeqcpp/PhreeqcKeywords" -) - -enable_testing() - -add_test(NAME run_pytest - COMMAND ${PYTHON_EXECUTABLE} -m pytest - WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/tests -) - -set_tests_properties(run_pytest - PROPERTIES ENVIRONMENT "PYTHONPATH=$") diff --git a/pybind/docstrings.h b/pybind/docstrings.h deleted file mode 100644 index fbbe1e38a..000000000 --- a/pybind/docstrings.h +++ /dev/null @@ -1,124 +0,0 @@ -#ifndef _INC_DOCSTRINGS_H -#define _INC_DOCSTRINGS_H - -#include - -// get_component_name_docstring -// get_components_docstring -// get_grid_size_docstring -// initialize_docstring -// update_docstring -// update_until_docstring - -/////////////////////////////////////////////////////////////////////////////// - -const std::string get_component_name_docstring = -R"pbdoc(Name of the component - -get_component_name returns the component name--"BMI PhreeqcRM". -BMI PhreeqcRM is a partial interface to PhreeqcRM, and provides -the methods to implement chemical reactions in a multicomponent -transport model. All of the native PhreeqcRM methods (non BMI -methods) are also available, which provides a complete -interface to PhreeqcRM. - -Returns -------- -str - The string "BMI PhreeqcRM". - -Examples --------- - ->>> import phreeqcrm as rm ->>> bmi = rm.bmi_phreeqcrm() ->>> print(bmi.get_component_name()) -BMI PhreeqcRM -)pbdoc"; - -/////////////////////////////////////////////////////////////////////////////// - -const std::string get_components_docstring = -R"pbdoc(Component list that was generated by calls to @FindComponents@. - -Returns -------- -tuple(str) - Each str is a component name. - -Examples --------- - ->>> import phreeqcrm as rm ->>> bmi = rm.bmi_phreeqcrm() ->>> bmi.initialize("AdvectBMI_py.yaml") ->>> components = bmi.get_components() ->>> print(components) -('H', 'O', 'Charge', 'Ca', 'Cl', 'K', 'N', 'Na') -)pbdoc"; - -/////////////////////////////////////////////////////////////////////////////// - -const std::string get_grid_size_docstring = -R"pbdoc(Get the total number of cells defined. - -get_grid_size returns the number of cells specified -by the YAML configuration file. If not specified a -default value of 10 is used. - -Parameters ----------- -grid : int - Grid number, only grid 0 is considered. - -Returns -------- -int - Same value as GetGridCellCount is returned for grid 0; - 0 for all other values of `grid`. -)pbdoc"; - -/////////////////////////////////////////////////////////////////////////////// - -const std::string initialize_docstring = -R"pbdoc(Initializes the bmi_phreeqcrm instance. - -A YAML file used for initialization contains a YAML map of -bmi_phreeqcrm methods and the arguments corresponding to -each method. - -Parameters ----------- -config_file : str, optional - The path to the model configuration file. -)pbdoc"; - -/////////////////////////////////////////////////////////////////////////////// - -const std::string update_docstring = -R"pbdoc(Advance model state by one time step. - -Runs PhreeqcRM for one time step. This method is equivalent to -@ref RunCells. PhreeqcRM will equilibrate the solutions with all equilibrium -reactants (EQUILIBRIUM_PHASES, EXCHANGE, GAS_PHASE, SOLID_SOLUTIONS, and SURFACE) -and integrate KINETICS reactions for the specified time step -(@ref SetValue "TimeStep" or @ref SetTimeStep). -)pbdoc"; - -/////////////////////////////////////////////////////////////////////////////// - -const std::string update_until_docstring = -R"pbdoc(Advance model state until the given time. - -Parameters ----------- -end_time : float - Time at the end of the time step. - -update_until is the same as @ref update, except the time step is calculated -from the argument @a end_time. The time step is calculated to be @a end_time minus -the current time (@ref GetCurrentTime). -)pbdoc"; - - -#endif /* _INC_DOCSTRINGS_H */ \ No newline at end of file diff --git a/pybind/pybind.cpp b/pybind/pybind.cpp deleted file mode 100644 index 590f7a0e0..000000000 --- a/pybind/pybind.cpp +++ /dev/null @@ -1,2211 +0,0 @@ -#include -#include -#include - -#include -#include -#include -// #include -// #include - -#include "BMIPhreeqcRM.h" -#include "VarManager.h" - -#define STRINGIFY(x) #x -#define MACRO_STRINGIFY(x) STRINGIFY(x) - -namespace py = pybind11; - - - -//{{ -#define PYBIND11_DETAILED_ERROR_MESSAGES -//}} - -//#include -//#include -//#include -//#include -//{{ -//#include -//#include -//#include -//#include -//#include -//}} - -namespace py = pybind11; - -// Function that returns a shared_ptr to a std::vector with some values -std::shared_ptr> get_shared_vector() { - auto v = std::make_shared>(std::initializer_list{1.0, 2.0, 3.0, 4.0, 5.0}); - return v; -} - -// Function that wraps a numpy array around a shared_ptr to a std::vector -py::array_t wrap_shared_vector(const std::shared_ptr>& v) { - auto size = v->size(); - auto data = v->data(); - py::capsule free_when_done(data, [](void* f) { - double* data = reinterpret_cast(f); - std::shared_ptr>* v = reinterpret_cast>*>(data); - delete v; - }); - return py::array(size, data, free_when_done); -} - -//// Module definition -//PYBIND11_MODULE(phreeqcrm, m) { -// //m.def("get_shared_vector", &get_shared_vector, "Get a shared_ptr to a std::vector"); -// //m.def("wrap_shared_vector", &wrap_shared_vector, "Wrap a numpy array around a shared_ptr to a std::vector"); -// m.def("echo3", [](py::str s) { -// if (s.is_none()) std::cout << "empty"; -// else std::cout << s; -// }, "A function that outputs a string.", py::arg("s").none(true) = py::none()); -//} - -//template -//py::array_t get_value(std::string name, py::array_t array) -//{ -// auto r = array.mutable_unchecked<1>(); -// for (py::ssize_t i = 0; i < r.shape(0); i++) { -// r(i) = 2 * r(i); -// } -// return array; -//} - - -py::array_t use_mutable_unchecked(std::string name, py::array_t array) -{ - auto r = array.mutable_unchecked<1>(); - for (py::ssize_t i = 0; i < r.shape(0); i++) { - r(i) = 2 * r(i); - } - return array; -} - -py::array get_value(std::string name, py::array array) -{ - py::buffer_info buffer = array.request(); - std::cout << "itemsize = " << buffer.itemsize << std::endl; - std::cout << "size = " << buffer.size << std::endl; - std::cout << "format = " << buffer.format << std::endl; - std::cout << "ndim = " << buffer.ndim << std::endl; - std::cout << "shape.size() = " << buffer.shape.size() << std::endl; - - std::cout << "shape: ("; - for (auto dim : buffer.shape) { - std::cout << dim << ", "; - } - std::cout << ")" << std::endl; - - - std::cout << "strides.size() = " << buffer.strides.size() << std::endl; - std::cout << "readonly = " << buffer.readonly << std::endl; - - if (buffer.ndim != 1) { - throw std::domain_error("array has incorrect number of dimensions: " - + std::to_string(buffer.ndim) + "; expected " - + std::to_string(1)); - } - - return array; -} - -py::array get_value_ptr(std::string name) -{ - // numpy.float64 - static std::vector v = { 1.2, 2.2, 3.2, 4.2, 5.2 }; - - py::buffer_info buf( - v.data(), - sizeof(double), - py::format_descriptor::format(), - 1, - { v.size() }, - { sizeof(double) } - ); - - ///py::detail::get_numpy_internals().get_type_info(true)->format_str; - - ///py::dtype numpy_dtype_double(py::format_descriptor::format()); - // numpy_dtype_double.type().str() - - std::cout << py::format_descriptor::format() << " = " << py::detail::get_numpy_internals().get_type_info(true)->format_str << std::endl; - std::cout << py::format_descriptor::format() << " = " << py::detail::get_numpy_internals().get_type_info(true)->format_str << std::endl; - - return py::array_t(buf); -} - -//std::string get_dtype_string(py::array_t& arr) { -// py::dtype dtype = arr.dtype(); // Get the dtype of the input array -// py::dtype::type type = dtype.type(); // Get the numpy dtype object for the dtype -// return type.str(); // Get the numpy dtype string for the numpy dtype object -//} - -// this requires #include -//std::string double_to_numpy_dtype() { -// std::string c_format = py::format_descriptor::format(); // Get the C type format string for double -// PyObject* numpy_dtype = PyArray_DescrFromType(NPY_NOTYPE); // Create a NumPy dtype object -// PyObject* resolved_dtype = PyArray_DescrFromTypestr(c_format.c_str()); // Resolve the NumPy dtype from the C type format string -// PyObject* dtype_repr = PyObject_Repr(resolved_dtype); // Get the string representation of the resolved dtype -// std::string numpy_dtype_str = PyUnicode_AsUTF8(dtype_repr); // Convert the dtype representation to a C++ string -// Py_XDECREF(numpy_dtype); -// Py_XDECREF(resolved_dtype); -// Py_XDECREF(dtype_repr); -// return numpy_dtype_str; -//} - -//py::array BMIPhreeqcRM::get_value_ptr(std::string name) -//{ -// // numpy.float64 -// //static std::vector v = { 1.2, 2.2, 3.2, 4.2, 5.2 }; -// -// //assert(name == "Temperature"); -// //std::vector v; -// //this->GetValue(name, v); -// -// -// RMVARS v_enum = this->var_man->GetEnum(name); -// if (v_enum == RMVARS::NotFound) throw std::runtime_error("Unknown variable name"); -// -// BMIVariant& bv = this->var_man->VariantMap[v_enum]; -// -// // int dim = bv.GetDim(); -// // std::cout << "0 dim = " << dim << std::endl; -// -// bool read_only = !bv.GetHasSetter(); -// -// //assert(name == "ComponentCount"); -// //std::vector v; -// //// this->GetValue(name, v); -// //int n; -// //this->GetValue(name, n); -// //v.push_back(n); -// -// void* ptr = this->GetValuePtr(name); // GetValuePtr takes a variable name and returns a pointer to a current copy of the variable values.Unlike the buffer returned from @ref GetValue, the reference always points to the current values of the variable, even if the model's state has changed. -// int nbytes = this->GetVarNbytes(name); // retrieves the total number of bytes that are set for a variable with @ref SetValue -// int itemsize = this->GetVarItemsize(name); // GetVarItemsize retrieves size of an individual item that can be set or retrived.Sizes may be sizeof(int), sizeof(double), or a character length for string variables. -// std::string type = this->GetVarType(name); // std::string GetVarType(const std::string name) -// -// // betting that GetDim call needs to be after first call to GetValuePtr() -// int dim = bv.GetDim(); -// std::cout << "0 dim = " << dim << std::endl; -// -// -// //ssize_t itemsize{ sizeof(double) }; -// //std::string format(py::format_descriptor::format()); -// std::string format; -// std::vector strides; -// if (type == "float64") -// { -// assert(itemsize == sizeof(double)); -// format = py::format_descriptor::format(); -// strides.push_back(sizeof(double)); -// } -// else if (type == "int32") -// { -// assert(itemsize == sizeof(int)); -// format = py::format_descriptor::format(); -// strides.push_back(sizeof(int)); -// } -// else -// { -// assert(false); -// throw std::runtime_error("Unknown dtype"); -// } -// -// int count = nbytes / itemsize; -// //assert(count == v.size()); -// std::cout << "count = " << count << std::endl; -// std::cout << "dim = " << dim << std::endl; -// assert(count == dim); -// -// //buffer_info_format; -// -// int nn = *((int*)ptr); -// -// //py::buffer_info buf( -// // ptr, -// // itemsize, // sizeof(double), -// // format, // py::format_descriptor::format(), -// // 1, -// // { count }, -// // strides, //{ sizeof(int) } // { sizeof(double) } -// // read_only -// //); -// -// py::buffer_info buf( -// ptr, -// itemsize, // sizeof(double), -// format, // py::format_descriptor::format(), -// 1, -// { count }, -// strides //{ sizeof(int) } // { sizeof(double) } -// ); -// -// -// ///py::detail::get_numpy_internals().get_type_info(true)->format_str; -// -// ///py::dtype numpy_dtype_double(py::format_descriptor::format()); -// // numpy_dtype_double.type().str() -// -// //std::cout << py::format_descriptor::format() << " = " << py::detail::get_numpy_internals().get_type_info(true)->format_str << std::endl; -// //std::cout << py::format_descriptor::format() << " = " << py::detail::get_numpy_internals().get_type_info(true)->format_str << std::endl; -// -// if (type == "int") -// { -// //if (read_only) return py::array_t(buf, py::capsule(ptr, [](void* ptr) { /* no-op */})); -// -// //if (read_only) { -// // py::array_t a = py::array_t(buf); -// -// // bool writeable = a.writeable(); -// -// // //reinterpret_cast(a.ptr())->flags &= ~pybind11::detail::npy_api::NPY_ARRAY_WRITEABLE_; -// // reinterpret_cast(a.ptr())->flags &= ~py::detail::npy_api::NPY_ARRAY_WRITEABLE_; -// -// // writeable = a.writeable(); -// -// // return a; -// //} -// -// -// if (read_only) { -// py::array_t a = py::array_t(buf); -// assert(a.writeable()); -// py::detail::array_proxy(a.ptr())->flags &= ~py::detail::npy_api::NPY_ARRAY_WRITEABLE_; -// assert(!a.writeable()); -// //reinterpret_cast(a.ptr())->flags &= ~py::detail::npy_api::NPY_ARRAY_WRITEABLE_; -// return a; -// } -// -// return py::array_t(buf); -// } -// -// if (read_only) { -// py::array_t a = py::array_t(buf); -// reinterpret_cast(a.ptr())->flags &= ~py::detail::npy_api::NPY_ARRAY_WRITEABLE_; -// return a; -// } -// return py::array_t(buf); -//} - -//py::array BMIPhreeqcRM::get_value_ptr(std::string name) -//{ -// RMVARS v_enum = this->var_man->GetEnum(name); -// if (v_enum == RMVARS::NotFound) throw std::runtime_error("Unknown variable name"); -// -// BMIVariant& bv = this->var_man->VariantMap[v_enum]; -// // if (!bv.GetHasPtr()) throw std::runtime_error("This variable does not support get_value_ptr."); // this throws on "ComponentCount" -// -// -// void* ptr = this->GetValuePtr(name); // GetValuePtr takes a variable name and returns a pointer to a current copy of the variable values.Unlike the buffer returned from @ref GetValue, the reference always points to the current values of the variable, even if the model's state has changed. -// int nbytes = this->GetVarNbytes(name); // retrieves the total number of bytes that are set for a variable with @ref SetValue -// int itemsize = this->GetVarItemsize(name); // GetVarItemsize retrieves size of an individual item that can be set or retrived.Sizes may be sizeof(int), sizeof(double), or a character length for string variables. -// std::string type = this->GetVarType(name); // std::string GetVarType(const std::string name) -// -// // must call GetDim after GetValuePtr -// int dim = bv.GetDim(); -// -// bool read_only = !bv.GetHasSetter(); -// -// std::string format; -// std::vector strides; -// if (type == "float64") -// { -// assert(itemsize == sizeof(double)); -// format = py::format_descriptor::format(); -// strides.push_back(sizeof(double)); -// } -// else if (type == "int32") -// { -// assert(itemsize == sizeof(int)); -// format = py::format_descriptor::format(); -// strides.push_back(sizeof(int)); -// } -// else -// { -// assert(false); -// throw std::runtime_error("Unknown dtype"); -// } -// -// int count = nbytes / itemsize; -// assert(count == dim); -// -// py::buffer_info buf( -// ptr, -// itemsize, -// format, -// 1, -// { count }, -// strides -// ); -// -// -// if (type == "int32") -// { -// if (read_only) { -// py::array_t a = py::array_t(buf); -// assert(a.writeable()); -// py::detail::array_proxy(a.ptr())->flags &= ~py::detail::npy_api::NPY_ARRAY_WRITEABLE_; -// assert(!a.writeable()); -// return a; -// } -// return py::array_t(buf); -// } -// -// if (read_only) { -// py::array_t a = py::array_t(buf); -// assert(a.writeable()); -// py::detail::array_proxy(a.ptr())->flags &= ~py::detail::npy_api::NPY_ARRAY_WRITEABLE_; -// assert(!a.writeable()); -// return a; -// } -// -// -// //return py::array_t(buf); -// py::array_t a = py::array_t(buf); -// -// //auto r = a.mutable_unchecked<1>(); -// ////r.data(0); -// //assert(r.data(0) == ptr); -// //for (py::ssize_t i = 0; i < r.shape(0); i++) { -// // r(i) = 2 * r(i); -// //} -// -// -// //if (true) { -// // py::array_t a = py::array_t(buf); -// // assert(a.owndata()); -// // py::detail::array_proxy(a.ptr())->flags &= ~py::detail::npy_api::NPY_ARRAY_OWNDATA_; -// // assert(!a.owndata()); -// // return a; -// // // auto view = py::array_t::ensure(arr.ptr(), nullptr); -// //} -// -// //if (true) { -// // // py::array_t::ensure() -// // py::array_t arr({ count }, (const double*)ptr, ); -// // assert(arr.data() == ptr); -// // return arr; -// //} -// -// if (true) { -// py::array_t arr({ 2, 3 }); -// auto view = py::array_t::ensure(arr.ptr()); -// // assert(!view.owndata()); -// assert(arr.data() == view.data()); -// } -// -// // std::vector vec = { 1.2, 1.3, 1.4 }; -// if (type == "float64") { -// //// Get a pointer to the vector data. -// //const double* data_ptr = vec.data(); -// -// //// Get the size of the vector. -// ////const size_t size = vec.size(); -// //const ssize_t size = vec.size(); -// -// //// Create a NumPy array that shares its data with the vector. -// //py::array_t arr({ size }, data_ptr); -// -// // Check that the array shares its data with the vector. -// //assert(arr.data() == data_ptr); -// -// std::vector &vec = bv.GetDoubleVectorRef(); -// -// ///py::cast(vec, pybind11::return_value_policy::reference); -// -// //py::array_t arr( { static_cast(vec.size()) }, vec.data(), py::cast(vec, pybind11::return_value_policy::reference)); -// -// //py::array_t arr({ static_cast(vec.size()) }, vec.data(), py::cast(ptr, pybind11::return_value_policy::reference)); -// -// assert(ptr == bv.GetVoidPtr()); -// -// //py::array_t arr({ static_cast(vec.size()) }, vec.data(), py::cast(ptr, pybind11::return_value_policy::reference)); -// //py::array_t arr({ bv.GetDim() }, vec.data(), py::cast(ptr, pybind11::return_value_policy::reference)); -// ///py::array_t arr({ bv.GetDim() }, (double*)bv.GetVoidPtr(), py::cast(bv.GetVoidPtr(), pybind11::return_value_policy::reference)); -// py::array_t arr({ bv.GetDim() }, (double*)bv.GetVoidPtr(), py::cast(bv.GetVoidPtr())); -// -// assert(arr.data() == vec.data()); -// assert(!arr.owndata()); -// -// py::dtype dt = arr.dtype(); -// -// std::cout << "Array data type : " << dt.str() << std::endl; -// -// -// if (read_only) { -// assert(arr.writeable()); -// py::detail::array_proxy(arr.ptr())->flags &= ~py::detail::npy_api::NPY_ARRAY_WRITEABLE_; -// assert(!arr.writeable()); -// return arr; -// } -// -// return arr; -// } -// -// return a; -//} - - -//void copy_array(const py::array_t& src, const py::array_t& dst) { -// auto src_buf = src.request(); // obtain buffer information for source array -// auto dst_buf = dst.request(); // obtain buffer information for destination array -// double* src_ptr = static_cast(src_buf.ptr); // get pointer to source data -// double* dst_ptr = static_cast(dst.mutable_data()); // get pointer to mutable destination data -// std::memcpy(dst_ptr, src_ptr, src_buf.size * sizeof(double)); // copy data from source to destination -//} - -/* -* Get a copy of values of the given variable -*/ -py::array BMIPhreeqcRM::get_value(std::string name, py::array dest) -{ - if (!this->_initialized) throw NotIntialized(); - - assert(this->var_man); - RMVARS v_enum = this->var_man->GetEnum(name); - if (v_enum == RMVARS::NotFound) throw std::runtime_error("Unknown variable name"); - - BMIVariant& bv = this->var_man->VariantMap[v_enum]; - ///if (!bv.GetHasPtr()) throw std::runtime_error("This variable does not support get_value_ptr."); // this throws on "ComponentCount" - - assert(this->language == "Py"); - - // this call is REQUIRED for the BMIVariant::GetVoidPtr() to be valid - this->GetValuePtr(name); - - assert(this->GetVarType(name) == bv.GetPType()); - - // dest dtype - py::ssize_t dst_ndim = dest.ndim(); - py::dtype dst_dt = dest.dtype(); - std::string dst_dt_str = std::string(dst_dt.str()); - - // src array - auto src = get_value_ptr(name); - - // src dtype - py::ssize_t src_ndim = src.ndim(); - py::dtype src_dt = src.dtype(); - std::string src_dt_str = std::string(src_dt.str()); - - if (dst_dt != src_dt) - { - throw std::runtime_error("bad array dtype"); - } - - auto src_buf = src.request(); // obtain buffer information for source array - auto dst_buf = dest.request(); // obtain buffer information for destination array - void* src_ptr = src_buf.ptr; - void* dst_ptr = dest.mutable_data(); - - std::memcpy(dst_ptr, src_ptr, src_buf.size * src_buf.itemsize); - - return dest; -} - -py::array BMIPhreeqcRM::get_value_ptr(std::string name) -{ - if (!this->_initialized) throw NotIntialized(); - - assert(this->var_man); - RMVARS v_enum = this->var_man->GetEnum(name); - if (v_enum == RMVARS::NotFound) throw std::runtime_error("Unknown variable name"); - - BMIVariant& bv = this->var_man->VariantMap[v_enum]; - - //{{ - if (bv.HasPyArray()) - { - return bv.GetPyArray(); - } - //}} - - ///if (!bv.GetHasPtr()) throw std::runtime_error("This variable does not support get_value_ptr."); // this throws on "ComponentCount" - - assert(this->language == "Py"); - - // this call is REQUIRED for the BMIVariant::GetVoidPtr() to be valid - void* ptr = this->GetValuePtr(name); - - assert(this->GetVarType(name) == bv.GetPType()); - - bool read_only = !(bv.GetHasSetter()); - - py::array a = py::none(); - - if (bv.GetPType() == "float64") - { - assert(bv.GetItemsize() == sizeof(double)); - assert(bv.GetVoidPtr() == this->GetValuePtr(name)); - - py::array_t arr({ bv.GetDim() }, (double*)bv.GetVoidPtr(), py::cast(bv.GetVoidPtr())); - a = arr; - - // Check that the array shares its data - assert(!arr.owndata()); - assert(arr.data() == bv.GetVoidPtr()); - - py::dtype dt = arr.dtype(); - assert(std::string(dt.str()) == bv.GetPType()); - } - else if (bv.GetPType() == "int32") - { - assert(bv.GetItemsize() == sizeof(int)); - assert(bv.GetVoidPtr() == this->GetValuePtr(name)); - - py::array_t arr({ bv.GetDim() }, (int*)bv.GetVoidPtr(), py::cast(bv.GetVoidPtr())); - a = arr; - - // Check that the array shares its data - assert(!arr.owndata()); - assert(arr.data() == bv.GetVoidPtr()); - - py::dtype dt = arr.dtype(); - assert(std::string(dt.str()) == bv.GetPType()); - } - else if (bv.GetPType().rfind(" &strings = bv.GetStringVectorRef(); - assert(strings.size() > 0); - - a = py::cast(strings); - // string arrays are always readonly - read_only = true; - - py::dtype dt = a.dtype(); - assert(std::string(dt.str()) == bv.GetPType()); - } - else - { - assert(false); - } - - if (read_only) - { - assert(a.writeable()); - py::detail::array_proxy(a.ptr())->flags &= ~py::detail::npy_api::NPY_ARRAY_WRITEABLE_; - assert(!a.writeable()); - } - return bv.SetPyArray(a); -} - -py::array BMIPhreeqcRM::get_value_at_indices(std::string name, py::array dest, py::array indices) -{ - // name and initilized validated in get_value_ptr - py::array src = get_value_ptr(name); - - // validate dest and indices - py::buffer_info dest_info = dest.request(); - py::buffer_info indices_info = indices.request(); - - if (dest_info.ndim != 1 || indices_info.ndim != 1) - { - throw std::runtime_error("Both dest and indices must be 1-dimensional arrays."); - } - - if (dest_info.shape[0] != indices_info.shape[0]) - { - throw std::runtime_error("The size of dest and indices must be the same."); - } - - // src dtype - py::dtype src_dt = src.dtype(); - - // dest dtype - py::dtype dst_dt = dest.dtype(); - - // indices dtype - py::dtype indices_dt = indices.dtype(); - - // validate dtypes - if (dst_dt != src_dt || std::string(indices_dt.str()) != "int32") - { - throw std::runtime_error("bad array dtype"); - } - - - int* indices_data = static_cast(indices_info.ptr); - - std::string src_dt_str = std::string(src_dt.str()); - if (src_dt_str == "float64") - { - double* dest_data = static_cast(dest_info.ptr); - - py::array_t double_array(src); - auto src_data = double_array.unchecked<1>(); - - for (py::ssize_t i = 0; i < dest_info.shape[0]; ++i) - { - dest_data[i] = src_data(static_cast(indices_data[i])); - } - } - else if (src_dt_str == "int32") - { - int* dest_data = static_cast(dest_info.ptr); - - py::array_t int_array(src); - auto src_data = int_array.unchecked<1>(); - - for (py::ssize_t i = 0; i < dest_info.shape[0]; ++i) - { - dest_data[i] = src_data(static_cast(indices_data[i])); - } - } - else if (src_dt_str.rfind(" str_array(src); - //auto src_data = int_array.unchecked<1>(); - - //for (py::ssize_t i = 0; i < dest_info.shape[0]; ++i) - //{ - // dest_data[i] = src_data(static_cast(indices_data[i])); - //} - assert(false); - throw std::runtime_error("str are unsupported in get_value_at_indices"); - } - else - { - assert(false); - throw std::runtime_error("bad array dtype"); - } - - - return dest; -} - -void BMIPhreeqcRM::set_value(std::string name, py::array arr) -{ - if (!this->_initialized) throw NotIntialized(); - - RMVARS v_enum = this->var_man->GetEnum(name); - if (v_enum == RMVARS::NotFound) throw std::runtime_error("Unknown variable name"); - - BMIVariant& bv = this->var_man->VariantMap[v_enum]; - - if (!bv.GetInitialized()) - { - this->var_man->task = VarManager::VAR_TASKS::Info; - ((*this->var_man).*bv.GetFn())(); - } - - py::ssize_t ndim = arr.ndim(); - py::dtype dt = arr.dtype(); - std::string dt_str = std::string(dt.str()); - - std::string bv_str = bv.GetPType(); - int dim = bv.GetDim(); - - const py::ssize_t* shape = arr.shape(); - py::ssize_t sz = arr.size(); - - // Check dimension - if (!(ndim > 0 && shape[0] == dim)) - { - std::ostringstream oss; - oss << "dimension error in set_value: " << name; - this->ErrorMessage(oss.str()); - throw std::runtime_error(oss.str()); - } - - if (dt_str == "float64") - { - //std::vector& ref = this->var_man->VarExchange.GetDoubleVectorRef(); - //ref.resize(arr.size()); - //std::memcpy(ref.data(), arr.data(), arr.size() * sizeof(double)); - assert(bv.GetCType() == "double" && dim > 1); - this->SetValue(name, (void*)arr.data()); - } - else if (dt_str == "int32") - { - //std::vector& ref = this->var_man->VarExchange.GetIntVectorRef(); - //ref.resize(arr.size()); - //std::memcpy(ref.data(), arr.data(), arr.size() * sizeof(int)); - assert(bv.GetCType() == "int" && dim > 1); - this->SetValue(name, (void*)arr.data()); - } - else if (dt_str.rfind(" result; - for (auto item : arr) - { - result.push_back(py::cast(item)); - } - - size_t n = result.size(); - assert(n > 0); - assert(shape[0] == n); - assert(ndim == 1); - assert(sz == n); - - if (n == 1) - { - this->SetValue(name, result[0]); - - std::string value; - this->GetValue(name, value); - assert(value == "prefix"); - } - else - { - // this should have been caught in check dims above - assert(false); - } - } - else if (dt_str == "bool") - { - std::vector result; - for (auto item : arr) - { - bool b = py::cast(item); - result.push_back(py::cast(item)); - } - - size_t n = result.size(); - assert(n > 0); - assert(shape[0] == n); - assert(ndim == 1); - assert(sz == n); - - if (n == 1) - { - this->SetValue(name, result[0]); - } - else - { - // this should have been caught in check dims above - assert(false); - } - } - return; -} - - -void BMIPhreeqcRM::set_value_at_indices(std::string name, py::array indices, py::array src) -{ - // name and initilized validated in get_value_ptr - py::array dest = get_value_ptr(name); - //{{ - py::buffer_info dest_info = dest.request(); - //}} - - // validate src and indices - py::buffer_info src_info = src.request(); - py::buffer_info indices_info = indices.request(); - - if (src_info.ndim != 1 || indices_info.ndim != 1) - { - throw std::runtime_error("Both src and indices must be 1-dimensional arrays."); - } - - if (src_info.shape[0] != indices_info.shape[0]) - { - throw std::runtime_error("The size of src and indices must be the same."); - } - - // src dtype - py::dtype src_dt = src.dtype(); - - // dest dtype - py::dtype dst_dt = dest.dtype(); - - // indices dtype - py::dtype indices_dt = indices.dtype(); - - // validate dtypes - if (dst_dt != src_dt || std::string(indices_dt.str()) != "int32") - { - throw std::runtime_error("bad array dtype"); - } - - int* indices_data = static_cast(indices_info.ptr); - - std::string src_dt_str = std::string(src_dt.str()); - if (src_dt_str == "float64") - { - double* dest_data = static_cast(dest_info.ptr); - - py::array_t double_array(src); - auto src_data = double_array.unchecked<1>(); - - for (py::ssize_t i = 0; i < src_info.shape[0]; ++i) - { - dest_data[indices_data[i]] = src_data(static_cast(i)); - } - } - else if (src_dt_str == "int32") - { - int* dest_data = static_cast(dest_info.ptr); - - py::array_t int_array(src); - auto src_data = int_array.unchecked<1>(); - - for (py::ssize_t i = 0; i < src_info.shape[0]; ++i) - { - dest_data[indices_data[i]] = src_data(static_cast(i)); - } - } - else if (src_dt_str.rfind(" str_array(src); - //auto src_data = int_array.unchecked<1>(); - - //for (py::ssize_t i = 0; i < dest_info.shape[0]; ++i) - //{ - // dest_data[i] = src_data(static_cast(indices_data[i])); - //} - assert(false); - throw std::runtime_error("str are unsupported in get_value_at_indices"); - } - else - { - assert(false); - throw std::runtime_error("bad array dtype"); - } -} - - -// def get_grid_shape(self, grid: int, shape : np.ndarray)->np.ndarray: - - -// see https://stackoverflow.com/questions/3898572/what-are-the-most-common-python-docstring-formats -// see https://www.sphinx-doc.org/en/master/usage/extensions/example_google.html#example-google -// see https://numpydoc.readthedocs.io/en/latest/example.html - - -// def foo(var1, var2, *args, long_var_name = "hi", only_seldom_used_keyword = 0, **kwargs) : -const std::string numpydoc_example_docstring = -R"pbdoc(Summarize the function in one line. - -Several sentences providing an extended description.Refer to -variables using back - ticks, e.g. `var`. - -Parameters ----------- -var1 : array_like - Array_like means all those objects -- lists, nested lists, etc. -- - that can be converted to an array.We can also refer to - variables like `var1`. -var2 : int - The type above can either refer to an actual Python type - (e.g. ``int``), or describe the type of the variable in more - detail, e.g. ``(N, ) ndarray`` or ``array_like``. -* args : iterable - Other arguments. -long_var_name : {'hi', 'ho'}, optional - Choices in brackets, default first when optional. - -Returns -------- -type - Explanation of anonymous return value of type ``type``. -describe : type - Explanation of return value named `describe`. -out : type - Explanation of `out`. -type_without_description - -Other Parameters ----------------- -only_seldom_used_keyword : int, optional - Infrequently used parameters can be described under this optional - section to prevent cluttering the Parameters section. -** kwargs : dict - Other infrequently used keyword arguments.Note that all keyword - arguments appearing after the first parameter specified under the - Other Parameters section, should also be described under this - section. - -Raises ------- -BadException - Because you shouldn't have done that. - -See Also --------- -numpy.array : Relationship(optional). -numpy.ndarray : Relationship(optional), which could be fairly long, in - which case the line wraps here. -numpy.dot, numpy.linalg.norm, numpy.eye - -Notes ------ -Notes about the implementation algorithm(if needed). - -This can have multiple paragraphs. - -You may include some math : - -..math::X(e^ { j\omega }) = x(n)e^ { -j\omega n } - -And even use a Greek symbol like : math:`\omega` inline. - -References ----------- -Cite the relevant literature, e.g.[1]_.You may also cite these -references in the notes section above. - -.. [1] O.McNoleg, "The integration of GIS, remote sensing, - expert systems and adaptive co - kriging for environmental habitat - modelling of the Highland Haggis using object - oriented, fuzzy - logic - and neural - network techniques, " Computers & Geosciences, vol. 22, - pp. 585 - 588, 1996. - -Examples --------- -These are written in doctest format, and should illustrate how to -use the function. - ->>> a = [1, 2, 3] ->>> print([x + 3 for x in a]) -[4, 5, 6] ->>> print("a\nb") -a -b -)pbdoc"; - - -const std::string Numpydoc_docstring = -R"pbdoc(My numpydoc description of a kind -of very exhautive numpydoc format docstring. - -Parameters ----------- -first : array_like - the 1st param name `first` -second : - the 2nd param - third : {'value', 'other'}, optional - the 3rd param, by default 'value' - -Returns -------- -string - a value in a string - -Raises ------- -KeyError - when a key error -OtherError - when an other error -)pbdoc"; - - -// Finalize closes any files open in the BMIPhreeqcRM instance. - -const std::string finalize_docstring = -R"pbdoc(Closes any files open in the bmi_phreeqcrm instance. -)pbdoc"; - - -const std::string initialize_docstring = -R"pbdoc(Initializes the bmi_phreeqcrm instance. - -A YAML file used for initialization contains a YAML map of -bmi_phreeqcrm methods and the arguments corresponding to -each method. - -Parameters ----------- -config_file : str, optional - The path to the model configuration file. -)pbdoc"; - - -const std::string get_components_docstring = -R"pbdoc(Component list that was generated by calls to @FindComponents@. - -Returns -------- -tuple(str) - Each str is a component name. - -Examples --------- - ->>> import phreeqcrm as rm ->>> bmi = rm.bmi_phreeqcrm() ->>> bmi.initialize("AdvectBMI_py.yaml") ->>> components = bmi.get_components() ->>> print(components) -('H', 'O', 'Charge', 'Ca', 'Cl', 'K', 'N', 'Na') -)pbdoc"; - - -const std::string get_component_name_docstring = -R"pbdoc(Name of the component - -get_component_name returns the component name--"BMI PhreeqcRM". -BMI PhreeqcRM is a partial interface to PhreeqcRM, and provides -the methods to implement chemical reactions in a multicomponent -transport model. All of the native PhreeqcRM methods (non BMI -methods) are also available, which provides a complete -interface to PhreeqcRM. - -Returns -------- -str - The string "BMI PhreeqcRM". - -Examples --------- - ->>> import phreeqcrm as rm ->>> bmi = rm.bmi_phreeqcrm() ->>> print(bmi.get_component_name()) -BMI PhreeqcRM -)pbdoc"; - -const std::string get_grid_size_docstring = -R"pbdoc(Get the total number of cells defined. - -get_grid_size returns the number of cells specified -by the YAML configuration file. If not specified a -default value of 10 is used. - -Parameters ----------- -grid : int - Grid number, only grid 0 is considered. - -Returns -------- -int - Same value as GetGridCellCount is returned for grid 0; - 0 for all other values of `grid`. -)pbdoc"; - -const std::string update_docstring = -R"pbdoc(Advance model state by one time step. - -Runs PhreeqcRM for one time step. This method is equivalent to -@ref RunCells. PhreeqcRM will equilibrate the solutions with all equilibrium -reactants (EQUILIBRIUM_PHASES, EXCHANGE, GAS_PHASE, SOLID_SOLUTIONS, and SURFACE) -and integrate KINETICS reactions for the specified time step -(@ref SetValue "TimeStep" or @ref SetTimeStep). -)pbdoc"; - -const std::string update_until_docstring = -R"pbdoc(Advance model state until the given time. - -Parameters ----------- -end_time : float - Time at the end of the time step. - -update_until is the same as @ref update, except the time step is calculated -from the argument @a end_time. The time step is calculated to be @a end_time minus -the current time (@ref GetCurrentTime). -)pbdoc"; - - -#include -#include - -namespace py = pybind11; - -// Define a templated function that accepts an array-like variable of ints or doubles -template -void my_function(py::array_t input_array) { - // Get a pointer to the data buffer - auto ptr = static_cast(input_array.request().ptr); - - // Get the shape of the array - auto shape = input_array.shape(); - - // Get the number of dimensions - auto ndim = input_array.ndim(); - - // Print some information about the array - std::cout << "Array shape: ("; - for (int i = 0; i < ndim; i++) { - std::cout << shape[i]; - if (i < ndim - 1) { - std::cout << ", "; - } - } - std::cout << ")" << std::endl; - - // Access the elements of the array - for (int i = 0; i < shape[0]; i++) { - for (int j = 0; j < shape[1]; j++) { - auto value = *(ptr + i * shape[1] + j); - std::cout << value << " "; - } - std::cout << std::endl; - } -} - -//// Define the Python module -//PYBIND11_MODULE(example, m) { -// m.def("my_function", &my_function, "A function that accepts an array-like variable of doubles"); -// m.def("my_function", &my_function, "A function that accepts an array-like variable of ints"); -//} - -//void my_func(py::array_t arr) { -void my_func(py::list strings) -{ - std::vector result; - for (auto item : strings) - { - std::string s = py::cast(item); - result.push_back(py::cast(item)); - } - //return result; -} - -py::sequence BMIPhreeqcRM::process_sequence(py::sequence seq) -{ - // Get the length of the sequence - ssize_t size = py::len(seq); - - py::handle type = seq.get_type(); - std::string string = std::string(type.str()); - - py::sequence return_val = seq; - - py::dtype int_dt = py::dtype::of(); - py::dtype double_dt = py::dtype::of(); - //py::dtype string_dt = py::dtype::of(); - - - if (py::isinstance(seq)) - { - py::list list(size); - //return_val = list; - - for (ssize_t i = 0; i < size; ++i) - { - // Access the elements of the sequence - py::object element = seq[i]; - - int modified_element = py::cast(element) * 2; - - //list.append(modified_element); - list[i] = modified_element; - } - return list; - } - else if (py::isinstance(seq)) - { - //py::tuple tuple(size); - py::tuple tuple(0); - return_val = tuple; - - for (ssize_t i = 0; i < size; ++i) - { - // Access the elements of the sequence - py::object element = seq[i]; - - int modified_element = py::cast(element) * 2; - - //tuple[i] = modified_element; - tuple = tuple + py::make_tuple(modified_element); - } - return tuple; - } - else if (py::isinstance(seq)) - { - //py::array a = seq.cast(py::array)(); - py::array arr = seq.cast(); - py::dtype dt = arr.dtype(); - std::string dt_str = std::string(dt.str()); - - py::object item = seq[0]; - - py::handle type = item.get_type(); - std::string string = std::string(type.str()); - - int int_element = py::cast(item); - - double double_element = py::cast(item); - - //if (py::isinstance(item)) - if (dt_str == "float64") - { - py::array_t d_array(size); - return_val = d_array; - } - //else if (py::isinstance(item)) - else if (dt_str == "int32") - { - py::array_t i_array(size); - return_val = i_array; - } - } - else - { - // unsupported sequence - return py::none(); - } - - //py::list mult_by_2; - - // Iterate over the sequence - //for (ssize_t i = 0; i < size; ++i) - //{ - // // Access the elements of the sequence - // py::object element = seq[i]; - - // int modified_element = py::cast(element) * 2; - - // mult_by_2.append(modified_element); - - // // Process the element - // // ... - // py::handle type = element.get_type(); - // std::string string = std::string(type.str()); - //} - //return mult_by_2; - return return_val; -} - -PYBIND11_MODULE(phreeqcrm, m) { - - m.def("my_function", &my_function, "A function that accepts an array-like variable of doubles"); - m.def("my_function", &my_function, "A function that accepts an array-like variable of ints"); - - m.def("my_func", - &my_func, - "A function that takes a NumPy array of strings." - ); - - py::class_(m, "bmi_phreeqcrm" , py::dynamic_attr()) - .def(py::init()) - - // - // testing - // - - .def("process_sequence", - &BMIPhreeqcRM::process_sequence, - py::arg("seq") - ) - - // - // Model control functions. - // - - // virtual void Initialize(std::string config_file) = 0; - // def initialize(self, config_file: str) -> None: - .def("initialize", - &BMIPhreeqcRM::Initialize, - initialize_docstring.c_str(), - py::arg("config_file") = py::str("") - ) - - // virtual void Update() = 0; - // def update(self) -> None: - .def("update", - [](BMIPhreeqcRM& self) { - if (!self._initialized) throw NotIntialized(); - self.Update(); - }, - update_docstring.c_str() - ) - - // virtual void UpdateUntil(double time) = 0; - // def update_until(self, time: float) -> None: - .def("update_until", - [](BMIPhreeqcRM& self, double time) { - if (!self._initialized) throw NotIntialized(); - self.UpdateUntil(time); - }, - update_until_docstring.c_str(), - py::arg("time") - ) - - // virtual void Finalize() = 0; - // def finalize(self) -> None: - .def("finalize", - [](BMIPhreeqcRM& self) { - self.Finalize(); - } - ) - - - // - // Model information functions. - // - - // virtual std::string GetComponentName() = 0; - // def get_component_name(self) -> str: - .def("get_component_name", - [](BMIPhreeqcRM& self) { - return self.GetComponentName(); - }, - get_component_name_docstring.c_str() - ) - - // virtual int GetInputItemCount() = 0; - // def get_input_item_count(self) -> int: - .def("get_input_item_count", - &BMIPhreeqcRM::GetInputItemCount - ) - - // virtual int GetOutputItemCount() = 0; - // def get_input_item_count(self) -> int: - .def("get_output_item_count", - &BMIPhreeqcRM::GetOutputItemCount - ) - - // virtual std::vector GetInputVarNames() = 0; - // def get_input_var_names(self) -> Tuple[str]: - .def("get_input_var_names", - [](BMIPhreeqcRM& self) { - if (!self._initialized) throw NotIntialized(); - auto output = self.GetInputVarNames(); - py::tuple out = py::cast(output); - return out; - } - ) - - // virtual std::vector GetOutputVarNames() = 0; - // def get_output_var_names(self) -> Tuple[str]: - .def("get_output_var_names", - [](BMIPhreeqcRM& self) { - if (!self._initialized) throw NotIntialized(); - auto output = self.GetOutputVarNames(); - py::tuple out = py::cast(output); - return out; - } - ) - - // - // Variable information functions - // - - // virtual int GetVarGrid(std::string name) = 0; - // def get_var_grid(self, name: str) -> int: - .def("get_var_grid", - &BMIPhreeqcRM::GetVarGrid, - py::arg("name") - ) - - // virtual std::string GetVarType(std::string name) = 0; - // def get_var_type(self, name: str) -> str: - .def("get_var_type", - &BMIPhreeqcRM::GetVarType, - py::arg("name") - ) - - // virtual std::string GetVarUnits(std::string name) = 0; - // def get_var_units(self, name: str) -> str: - .def("get_var_units", - &BMIPhreeqcRM::GetVarUnits, - py::arg("name") - ) - - // virtual int GetVarItemsize(std::string name) = 0; - - // virtual int GetVarNbytes(std::string name) = 0; - // def get_var_nbytes(self, name: str) -> int: - .def("get_var_nbytes", - &BMIPhreeqcRM::GetVarNbytes, - py::arg("name") - ) - - // virtual std::string GetVarLocation(std::string name) = 0; - // def get_var_location(self, name: str) -> str: - .def("get_var_location", - [](BMIPhreeqcRM& self, std::string name) { - if (!self._initialized) throw NotIntialized(); - return self.GetVarLocation(name); - }, - py::arg("name") - ) - - // virtual double GetCurrentTime() = 0; - // def get_current_time(self) -> float: - .def("get_current_time", - [](BMIPhreeqcRM& self) { - if (!self._initialized) throw NotIntialized(); - return self.GetCurrentTime(); - } - ) - - // virtual double GetStartTime() = 0; - // def get_start_time(self) -> float: - .def("get_start_time", - [](BMIPhreeqcRM& self) { - if (!self._initialized) throw NotIntialized(); - return self.GetStartTime(); - } - ) - - - // virtual double GetEndTime() = 0; - // def get_end_time(self) -> float: - .def("get_end_time", - [](BMIPhreeqcRM& self) { - if (!self._initialized) throw NotIntialized(); - return self.GetEndTime(); - } - ) - - - // virtual std::string GetTimeUnits() = 0; - // def get_time_units(self) -> str: - .def("get_time_units", - [](BMIPhreeqcRM& self) { - return self.GetTimeUnits(); - } - ) - - // virtual double GetTimeStep() = 0; - // def get_time_step(self) -> float: - .def("get_time_step", - [](BMIPhreeqcRM& self) { - return self.GetTimeStep(); - } - ) - - - // - // Variable getters - // - - // virtual void GetValue(std::string name, void *dest) = 0; - // def get_value(self, name: str, dest: np.ndarray) -> np.ndarray: - .def("get_value", - &BMIPhreeqcRM::get_value, - py::arg("name"), - py::arg("dest") - ) - - // virtual void *GetValuePtr(std::string name) = 0; - // def get_value_ptr(self, name: str) -> np.ndarray: - .def("get_value_ptr", - &BMIPhreeqcRM::get_value_ptr, - py::arg("name") - ) - - // virtual void GetValueAtIndices(std::string name, void *dest, int *inds, int count) = 0; - // def get_value_at_indices( - // self, name: str, dest : np.ndarray, inds : np.ndarray - // )->np.ndarray: - - //.def("get_value_at_indices", - // [](BMIPhreeqcRM& self, std::string name, py::array dest) { - // auto resultobj = self.get_value_ptr(name).attr("take"); - // return resultobj(); - // } - //) - .def("get_value_at_indices", - &BMIPhreeqcRM::get_value_at_indices, - py::arg("name"), - py::arg("dest"), - py::arg("inds") - ) - - - //.def("my_sqrt", - // [](BMIPhreeqcRM& self, double value) { - // auto math = py::module::import("math"); - // auto resultobj = math.attr("sqrt")(value); - // return resultobj.cast(); - // }, - // py::arg("value") - //) - - - // - // Variable setters - // - - // virtual void SetValue(std::string name, void *src) = 0; - // def set_value(self, name: str, src: np.ndarray) -> None: - .def("set_value", - &BMIPhreeqcRM::set_value, - py::arg("name"), - py::arg("src") - ) - - // virtual void SetValueAtIndices(std::string name, int *inds, int count, void *src) = 0; - // def set_value_at_indices( - // self, name: str, inds : np.ndarray, src : np.ndarray - // )->None: - .def("set_value_at_indices", - &BMIPhreeqcRM::set_value_at_indices, - py::arg("name"), - py::arg("inds"), - py::arg("src") - ) - - // - // Grid information functions - // - - // virtual int GetGridRank(const int grid) = 0; - // def get_grid_rank(self, grid: int) -> int: - .def("get_grid_rank", - &BMIPhreeqcRM::GetGridRank, - py::arg("grid") - ) - - // virtual int GetGridSize(const int grid) = 0; - // def get_grid_size(self, grid: int) -> int: - .def("get_grid_size", - [](BMIPhreeqcRM& self, int grid) { - if (!self._initialized) throw NotIntialized(); - return self.GetGridSize(grid); - }, - get_grid_size_docstring.c_str(), - py::arg("grid") - ) - - // virtual std::string GetGridType(const int grid) = 0; - // def get_grid_type(self, grid: int) -> str: - .def("get_grid_type", - &BMIPhreeqcRM::GetGridType, - py::arg("grid") - ) - - - - // void GetGridShape(const int grid, int* shape) override - // def get_grid_shape(self, grid: int, shape: np.ndarray) -> np.ndarray: - .def("get_grid_shape", - [](BMIPhreeqcRM& self, int grid, py::array_t shape) { - throw NotImplemented(); - if (!self._initialized) throw NotIntialized(); - - int* data = shape.mutable_data(); - py::array::ShapeContainer shapeContainer{1}; - shape.reshape(shapeContainer); - data[0] = self.GetGridCellCount(); - return shape; - } - , - //get_grid_shape_docstring.c_str(), - py::arg("grid"), - py::arg("shape") - ) - - // virtual void GetGridSpacing(const int grid, double *spacing) = 0; - // def get_grid_spacing(self, grid: int, spacing: np.ndarray) -> np.ndarray: - .def("get_grid_spacing", - [](BMIPhreeqcRM& self, int grid, py::array_t spacing) { - throw NotImplemented(); - if (!self._initialized) throw NotIntialized(); - return spacing; - } - , - //get_grid_spacing_docstring.c_str(), - py::arg("grid"), - py::arg("spacing") - ) - - // virtual void GetGridOrigin(const int grid, double *origin) = 0; - // def get_grid_origin(self, grid: int, origin: np.ndarray) -> np.ndarray: - .def("get_grid_origin", - [](BMIPhreeqcRM& self, int grid, py::array_t origin) { - throw NotImplemented(); - if (!self._initialized) throw NotIntialized(); - return origin; - } - , - //get_grid_origin_docstring.c_str(), - py::arg("grid"), - py::arg("origin") - ) - - // virtual void GetGridX(const int grid, double* x) = 0; - // def get_grid_x(self, grid: int, x: np.ndarray) -> np.ndarray: - .def("get_grid_x", - [](BMIPhreeqcRM& self, int grid, py::array x) { - throw NotImplemented(); - return x; - }, - //get_grid_x_docstring.c_str(), - py::arg("grid"), - py::arg("x") - ) - - // virtual void GetGridY(const int grid, double* y) = 0; - // def get_grid_y(self, grid: int, y: np.ndarray) -> np.ndarray: - .def("get_grid_y", - [](BMIPhreeqcRM& self, int grid, py::array y) { - throw NotImplemented(); - return y; - }, - //get_grid_y_docstring.c_str(), - py::arg("grid"), - py::arg("y") - ) - - // virtual void GetGridZ(const int grid, double* z) = 0; - // def get_grid_y(self, grid: int, z: np.ndarray) -> np.ndarray: - .def("get_grid_z", - [](BMIPhreeqcRM& self, int grid, py::array z) { - throw NotImplemented(); - return z; - }, - //get_grid_z_docstring.c_str(), - py::arg("grid"), - py::arg("z") - ) - - // virtual int GetGridNodeCount(const int grid) = 0; - // def get_grid_node_count(self, grid: int) -> int: - .def("get_grid_node_count", - [](BMIPhreeqcRM& self, int grid) { - //throw NotImplemented(); - return self.GetGridNodeCount(grid); - }, - // get_grid_node_count_docstring.c_str(), - py::arg("grid") - ) - - // virtual int GetGridEdgeCount(const int grid) = 0; - // def get_grid_edge_count(self, grid: int) -> int: - .def("get_grid_edge_count", - [](BMIPhreeqcRM& self, int grid) { - // throw NotImplemented(); - return self.GetGridEdgeCount(grid); - }, - // get_grid_edge_count_docstring.c_str(), - py::arg("grid") - ) - - // virtual int GetGridFaceCount(const int grid) = 0; - // def get_grid_face_count(self, grid: int) -> int: - .def("get_grid_face_count", - [](BMIPhreeqcRM& self, int grid) { - // throw NotImplemented(); - return self.GetGridFaceCount(grid); - }, - // get_grid_face_count_docstring.c_str(), - py::arg("grid") - ) - - - // virtual void GetGridEdgeNodes(const int grid, int* edge_nodes) = 0; - // def get_grid_edge_nodes(self, grid: int, edge_nodes: np.ndarray) -> np.ndarray: - .def("get_grid_edge_nodes", - [](BMIPhreeqcRM& self, int grid, py::array_t edge_nodes) { - throw NotImplemented(); - return edge_nodes; - }, - // get_grid_edge_nodes_docstring.c_str(), - py::arg("grid"), - py::arg("edge_nodes") - ) - - // virtual void GetGridFaceEdges(const int grid, int* face_edges) = 0; - // def get_grid_face_edges(self, grid: int, face_edges: np.ndarray) -> np.ndarray: - .def("get_grid_face_edges", - [](BMIPhreeqcRM& self, int grid, py::array_t face_edges) { - throw NotImplemented(); - return face_edges; - }, - // get_grid_face_edges_docstring.c_str(), - py::arg("grid"), - py::arg("face_edges") - ) - - // virtual void GetGridFaceNodes(const int grid, int* face_nodes) = 0; - // def get_grid_face_nodes(self, grid: int, face_nodes: np.ndarray) -> np.ndarray: - .def("get_grid_face_nodes", - [](BMIPhreeqcRM& self, int grid, py::array_t face_nodes) { - throw NotImplemented(); - return face_nodes; - }, - // get_grid_face_nodes_docstring.c_str(), - py::arg("grid"), - py::arg("face_nodes") - ) - - // virtual void GetGridNodesPerFace(const int grid, int* nodes_per_face) = 0; - // def get_grid_nodes_per_face( - // self, grid: int, nodes_per_face: np.ndarray - // ) -> np.ndarray: - .def("get_grid_nodes_per_face", - [](BMIPhreeqcRM& self, int grid, py::array_t nodes_per_face) { - throw NotImplemented(); - return nodes_per_face; - }, - // get_grid_nodes_per_face_docstring.c_str(), - py::arg("grid"), - py::arg("nodes_per_face") - ) - - // - // Extras - // - - .def("get_components", - [](BMIPhreeqcRM& self) { - // initialization not necessary - auto output = self.GetComponents(); - py::tuple out = py::cast(output); - return out; - }, - get_components_docstring.c_str() - ) - - - .def_readonly("_initialized", &BMIPhreeqcRM::_initialized) - - ; - - //m.def("get_value", &get_value, "Get a shared_ptr to a std::vector"); - //m.def("get_value", &get_value, "std::vector"); - //m.def("get_value", &get_value, "std::vector"); - //m.def("get_value", &get_value, ""); - //m.def("get_value_ptr", &get_value_ptr, ""); - //m.def("use_mutable_unchecked", &use_mutable_unchecked, ""); -} - - - - - -// PYBIND11_MODULE(phreeqcrm, m) { -// m.doc() = "I'm a docstring hehe"; -// m.def("some_fn_python_name", &some_fn); -// py::class_( -// m, "bmi_phreeqcrm" -// ) -// .def(py::init()) -// .def("initialize", &BMIPhreeqcRM::Initialize) -// ; -// } - -int add(int x, int y=0) { - return x + y; -} - -const std::string bmi_initialize_docstring = -R"pbdoc(A YAML file can be used to initialize an instance of bmi_phreeqcrm. - -Parameters ----------- -filename : str, optional - Path to name of input file. -)pbdoc"; - - -const std::string bmi_set_value_docstring = -R"pbdoc("Use the vector of concentrations (c) to set the moles of components in each -reaction cell. The volume of water in a cell is the product of porosity -(:meth:`SetPorosity`), saturation (:meth:`SetSaturation`), and reference volume -(:meth:`SetRepresentativeVolume`). The moles of each component are determined -by the volume of water and per liter concentrations. If concentration units -(:meth:`SetUnitsSolution`) are mass fraction, the density (as specified by -:meth:`SetDensity`) is used to convert from mass fraction to per mass per liter. - -:param c: a list or numpy array of doubles. - -Args: - c (double list, numpy.ndarray, or tuple): Component concentrations. Size of - vector is ncomps times nxyz, where ncomps is the number of components as - determined by :meth:`FindComponents` or :meth:`GetComponentCount` and nxyz - is the number of grid cells in the user model (:meth:`GetGridCellCount`). -)pbdoc"; - - -#include - -// https://github.com/Deltares/xmipy/blob/f3954fd888e42d24920adc58543185b52f80a6c7/xmipy/xmiwrapper.py#L313 - -//py::array BMIPhreeqcRM::get_value_test(std::string name, py::array dest/* = py::none()*/) -////py::array BMIPhreeqcRM::get_value_test(std::string name, py::array_t dest = py::none()) -//{ -// //{{ -// // make sure that optional array is of correct layout: -// if (!dest.is_none()) -// { -// if (dest.flags() & py::array::f_style) -// { -// throw std::runtime_error("Array should have C layout"); -// } -// } -// -// RMVARS v_enum = this->var_man->GetEnum(name); -// return py::none(); -// //return dest; -// -// -// //if (v_enum != RMVARS::NotFound) -// //{ -// // // first deal with scalars -// // BMIVariant& bv = this->var_man->VariantMap[v_enum]; -// -// // int dims = this->var_man->VarExchange.GetDim(); -// -// // if (dims > 0) -// // { -// // py::array out = py::cast(this->var_man->VarExchange.GetDoubleVectorRef()); -// // return out; -// // } -// // else -// // { -// // dest = this->var_man->VarExchange.GetDVar(); -// // double output = 3.14; -// // return py::cast(); -// // } -// -// //} -// -// -// ////}} -// -// -// //RMVARS v_enum = this->var_man->GetEnum(name); -// //if (v_enum != RMVARS::NotFound) -// //{ -// // BMIVariant& bv = this->var_man->VariantMap[v_enum]; -// // //VarManager::VarFunction fn = this->var_man->GetFn(v_enum); -// // if (!bv.GetInitialized()) -// // { -// // this->var_man->task = VarManager::VAR_TASKS::Info; -// // ((*this->var_man).*bv.GetFn())(); -// // } -// // this->var_man->task = VarManager::VAR_TASKS::GetVar; -// // ((*this->var_man).*bv.GetFn())(); -// // -// // int dims = this->var_man->VarExchange.GetDim(); -// -// // if (this->var_man->VarExchange.GetCType() == "double") -// // { -// // if (dims > 0) -// // { -// // py::array out = py::cast(this->var_man->VarExchange.GetDoubleVectorRef()); -// // return out; -// // } -// // else -// // { -// // dest = this->var_man->VarExchange.GetDVar(); -// // double output = 3.14; -// // return py::cast(); -// // } -// // } -// -// // //dest = this->var_man->VarExchange.GetDoubleVectorRef(); -// // //return; -// //} -// //return py::none(); -// -// -// ////if (name.compare("Pi") == 0) -// ////{ -// //// double output = 3.14; -// //// return py::cast(output); -// ////} -// ////else -// ////{ -// //// std::vector output(20, 3.14); -// //// py::array out = py::cast(output); -// //// return out; -// ////} -//} - -//#include -// -//void echo(std::string s) -//{ -// std::cout << s << std::endl; -//} -// -////int my_len(py::array_t a) -////{ -//// if (a.is_none()) return 0; -//// -//// py::buffer_info buf = a.request(); -//// return buf.shape[0]; -////} -// -//int my_len(py::array a) -//{ -// if (a.is_none()) return 0; -// -// py::buffer_info buf = a.request(); -// return buf.shape[0]; -//} -// -// -////using array_i = pybind11::array_t; -//using array_i = pybind11::array_t; -// -//struct myObj_array { -// array_i ids; -//}; -// - - - - -//PYBIND11_MODULE(phreeqcrm, m) { -// -// py::class_(m, "myObj_array_int") -// .def(py::init([](array_i py_ids) { -// std::cout << "ndim " << py_ids.ndim() << std::endl; -// if (py_ids.ndim()) -// return myObj_array{ py_ids }; -// else -// return myObj_array{ array_i(0) }; -// }), py::arg("py_ids") = (array_i*) nullptr -// ); -// -// -// -// //m.def("add", &add, "A function that adds two integers.", py::arg("x"), py::arg("y")=0); -// -// m.def("echo", &echo, "A function that outputs a string.", py::arg("s") = std::string("Hello")); -// -// m.def("echo2", [](py::str s) { -// if (s.is_none()) std::cout << "empty"; -// else std::cout << s; -// }, "A function that outputs a string.", py::arg("s") = static_cast(nullptr); -// -// -// m.def("echo3", [](py::object s) { -// if (s.is_none()) std::cout << "empty"; -// else std::cout << s; -// }, "A function that outputs a string.", py::arg("s") = py::none()); -// -// -// m.def("my_len", &my_len, "A function that outputs a string.", py::arg("a") = py::none()); -// -// m.def("my_len2", -// [](std::optional a) { -// if (!a.has_value()) { -// return 0; -// } -// return 10; -// }, -// py::arg("a") = py::none() -// ); -// -// -// -// m.doc() = R"pbdoc( -// PhreeqcRM plugin -// ----------------------- -// -// .. currentmodule:: phreeqcrm -// -// .. autosummary:: -// :toctree: _generate -// -// )pbdoc"; -// -// py::class_( -// m, "bmi_phreeqcrm" /* , py::dynamic_attr()*/ -// ) -// .def(py::init()) -// -// // Model control functions. -// // .def("initialize", &BMIPhreeqcRM::Initialize, R"pbdoc( -// // A YAML file can be used to initialize an instance of bmi_phreeqcrm. -// -// // Parameters -// // ---------- -// // filename : str, optional -// // Path to name of input file. -// // )pbdoc", py::arg("filename")) -// -// // .def("initialize2", [](BMIPhreeqcRM &self, py::str filename) { -// // if (filename.is_none()) { -// // return self.Initialize(""); -// // } -// // return self.Initialize(std::string(py::str(filename))); -// // }, bmi_initialize_docstring2.c_str(), py::arg("filename") = py::none()) -// -// .def("initialize", [](BMIPhreeqcRM &self, py::str filename) { -// if (filename.is_none()) { -// return self.Initialize(""); -// } -// return self.Initialize(std::string(py::str(filename))); -// }, bmi_initialize_docstring.c_str(), py::arg("filename") = py::none()) -// -// -// -// // Model information functions. -// .def("get_component_name", &BMIPhreeqcRM::GetComponentName) -// -// // .def("get_input_var_names", &BMIPhreeqcRM::GetInputVarNames) -// .def("get_input_var_names", [](BMIPhreeqcRM &self) { -// auto output = self.GetInputVarNames(); -// py::tuple out = py::cast(output); -// return out; -// }) -// -// .def("get_output_var_names", [](BMIPhreeqcRM &self) { -// auto output = self.GetOutputVarNames(); -// py::tuple out = py::cast(output); -// return out; -// }) -// -// -// -// // Variable getters -// // virtual void GetValue(std::string name, void *dest) = 0; -// .def("get_value", [](BMIPhreeqcRM &self, const std::string &input) { -// std::vector output; -// self.GetValue(input, output); -// py::array out = py::cast(output); -// return out; -// }) -// // virtual void *GetValuePtr(std::string name) = 0; -// // virtual void GetValueAtIndices(std::string name, void *dest, int *inds, int count) = 0; -// -// //.def("get_value_test", &BMIPhreeqcRM::get_value_test) -// //.def("get_value_test", &BMIPhreeqcRM::get_value_test, py::arg("name"), py::arg("dest") = py::none()) -// -// .def("get_value_test", [](BMIPhreeqcRM& self, const std::string& name, py::array dest) { -// std::vector output; -// //self.GetValue(input, output); -// //py::array out = py::cast(dest); -// //return out; -// return dest; -// }, py::arg("name"), py::arg("dest") = py::none()) -// -// -// -// // Variable setters -// // virtual void SetValue(std::string name, void *src) = 0; -// //.def("set_value", py::overload_cast>(&BMIPhreeqcRM::SetValue), py::arg("name", std::string()), py::arg("src", std::vector())) -// // .def("set_value", py::overload_cast>(&BMIPhreeqcRM::SetValue), py::arg("name"), py::arg("src").noconvert(), py::keep_alive<1, 2>(), -// // bmi_set_value_docstring.c_str() -// // ) -// -// -// .def("set_value", [](BMIPhreeqcRM &self, std::string name, const py::array_t& src) { -// // Extract the buffer information from the numpy array. -// py::buffer_info info = src.request(); -// -// // Construct a std::vector from the buffer data. -// double* ptr = static_cast(info.ptr); -// size_t size = info.shape[0]; -// std::vector src_vec(ptr, ptr + size); -// -// self.SetValue(name, src_vec); -// }, -// py::arg("name"), -// //py::arg("src").noconvert(), // this will fail if a List is used -// py::arg("src"), -// bmi_set_value_docstring.c_str() -// ) -// -// -// .def("set_value", [](BMIPhreeqcRM &self, std::string name, const py::array_t& src) { -// // Extract the buffer information from the numpy array. -// py::buffer_info info = src.request(); -// -// // Construct a std::vector from the buffer data. -// int* ptr = static_cast(info.ptr); -// size_t size = info.shape[0]; -// std::vector src_vec(ptr, ptr + size); -// -// self.SetValue(name, src_vec); -// }, -// py::arg("name"), -// //py::arg("src").noconvert(), // this will fail if a List is used -// py::arg("src"), -// bmi_set_value_docstring.c_str() -// ) -// -// -// -// -// // .def("set_value", [](BMIPhreeqcRM &self, const std::string &input, std::vector values) { -// // self.GetValue(input, output); -// // py::array out = py::cast(output); -// // return out; -// // }) -// -// -// // virtual void SetValueAtIndices(std::string name, int *inds, int count, void *src) = 0; -// -// -// // Grid information functions -// .def("get_grid_size", &BMIPhreeqcRM::GetGridSize) -// // virtual int GetGridRank(const int grid) = 0; -// // virtual int GetGridSize(const int grid) = 0; -// // virtual std::string GetGridType(const int grid) = 0; -// -// // virtual void GetGridShape(const int grid, int *shape) = 0; -// // virtual void GetGridSpacing(const int grid, double *spacing) = 0; -// // virtual void GetGridOrigin(const int grid, double *origin) = 0; -// -// // virtual void GetGridX(const int grid, double *x) = 0; -// // virtual void GetGridY(const int grid, double *y) = 0; -// // virtual void GetGridZ(const int grid, double *z) = 0; -// -// // virtual int GetGridNodeCount(const int grid) = 0; -// // virtual int GetGridEdgeCount(const int grid) = 0; -// // virtual int GetGridFaceCount(const int grid) = 0; -// -// // virtual void GetGridEdgeNodes(const int grid, int *edge_nodes) = 0; -// // virtual void GetGridFaceEdges(const int grid, int *face_edges) = 0; -// // virtual void GetGridFaceNodes(const int grid, int *face_nodes) = 0; -// // virtual void GetGridNodesPerFace(const int grid, int *nodes_per_face) = 0; -// -// -// // Extras -// .def("get_components", &BMIPhreeqcRM::GetComponents) -// -// ; -// -// -// -//#ifdef VERSION_INFO -// m.attr("__version__") = MACRO_STRINGIFY(VERSION_INFO); -//#else -// m.attr("__version__") = "dev"; -//#endif -//} - - - -// #include - -// #define STRINGIFY(x) #x -// #define MACRO_STRINGIFY(x) STRINGIFY(x) - -// int add(int i, int j) { -// return i + j; -// } - -// namespace py = pybind11; - -// PYBIND11_MODULE(cmake_example, m) { -// m.doc() = R"pbdoc( -// Pybind11 example plugin -// ----------------------- - -// .. currentmodule:: cmake_example - -// .. autosummary:: -// :toctree: _generate - -// add -// subtract -// )pbdoc"; - -// m.def("add", &add, R"pbdoc( -// Add two numbers - -// Some other explanation about the add function. -// )pbdoc"); - -// m.def("subtract", [](int i, int j) { return i - j; }, R"pbdoc( -// Subtract two numbers - -// Some other explanation about the subtract function. -// )pbdoc"); - -// #ifdef VERSION_INFO -// m.attr("__version__") = MACRO_STRINGIFY(VERSION_INFO); -// #else -// m.attr("__version__") = "dev"; -// #endif -// } diff --git a/pybind/pybind.cxx b/pybind/pybind.cxx deleted file mode 100644 index 578e1bffb..000000000 --- a/pybind/pybind.cxx +++ /dev/null @@ -1,986 +0,0 @@ -#include -#include -#include - -#include -#include -#include - -#include "BMIPhreeqcRM.h" -#include "VarManager.h" - -#include "docstrings.h" - -#define STRINGIFY(x) #x -#define MACRO_STRINGIFY(x) STRINGIFY(x) - -namespace py = pybind11; - -#define PYBIND11_DETAILED_ERROR_MESSAGES - -/* -* Get a copy of values of the given variable -*/ -py::array BMIPhreeqcRM::get_value(std::string name, py::array dest) -{ - if (!this->_initialized) throw NotIntialized(); - - assert(this->var_man); - RMVARS v_enum = this->var_man->GetEnum(name); - if (v_enum == RMVARS::NotFound) throw std::runtime_error("Unknown variable name"); - - BMIVariant& bv = this->var_man->VariantMap[v_enum]; - - assert(this->language == "Py"); - - if (bv.GetHasPtr()) - { - // this call is REQUIRED for the BMIVariant::GetVoidPtr() to be valid - this->GetValuePtr(name); - } - - assert(this->GetVarType(name) == bv.GetPType()); - - // dest dtype - py::ssize_t dst_ndim = dest.ndim(); - py::dtype dst_dt = dest.dtype(); - std::string dst_dt_str = std::string(dst_dt.str()); - - // check for scalars - int dim = bv.GetDim(); - if (!bv.GetHasPtr() && dim == 1) - { - // scalar - py::array a; - if (bv.GetPType() == "float64") - { - assert(bv.GetItemsize() == sizeof(double)); - assert(bv.GetNbytes() == sizeof(double)); - - void* dst_ptr = dest.mutable_data(); - const void* src_ptr = bv.GetDVarPtr(); - assert(src_ptr); - - std::memcpy(dst_ptr, src_ptr, bv.GetNbytes()); - - return dest; - } - else if (bv.GetPType() == "int32") - { - assert(bv.GetItemsize() == sizeof(int)); - assert(bv.GetNbytes() == sizeof(int)); - - void* dst_ptr = dest.mutable_data(); - const void* src_ptr = bv.GetIVarPtr(); - assert(src_ptr); - - std::memcpy(dst_ptr, src_ptr, bv.GetNbytes()); - - return dest; - } - else if (bv.GetPType().rfind(" strings; - strings.push_back(bv.GetStringVar()); - py::array src = py::cast(strings); - - ssize_t nsrc = src.request().itemsize; - ssize_t ndst = dest.request().itemsize; - - if (src.request().itemsize > dest.request().itemsize) - { - throw std::runtime_error("buffer too small"); - } - - std::memset(dest.mutable_data(), 0, dest.request().itemsize); - std::memcpy(dest.mutable_data(), src.request().ptr, src.request().itemsize); - - return dest; - } - else - { - assert(false); - } - } - - // src array - auto src = get_value_ptr(name); - - // src dtype - py::ssize_t src_ndim = src.ndim(); - py::dtype src_dt = src.dtype(); - std::string src_dt_str = std::string(src_dt.str()); - - if (dst_dt != src_dt) - { - throw std::runtime_error("bad array dtype"); - } - - auto src_buf = src.request(); // obtain buffer information for source array - auto dst_buf = dest.request(); // obtain buffer information for destination array - void* src_ptr = src_buf.ptr; - void* dst_ptr = dest.mutable_data(); - - std::memcpy(dst_ptr, src_ptr, src_buf.size * src_buf.itemsize); - - return dest; -} - -py::array BMIPhreeqcRM::get_value_at_indices(std::string name, py::array dest, py::array indices) -{ - // name and initilized validated in get_value_ptr - py::array src = get_value_ptr(name); - - // validate dest and indices - py::buffer_info dest_info = dest.request(); - py::buffer_info indices_info = indices.request(); - - if (dest_info.ndim != 1 || indices_info.ndim != 1) - { - throw std::runtime_error("Both dest and indices must be 1-dimensional arrays."); - } - - if (dest_info.shape[0] != indices_info.shape[0]) - { - throw std::runtime_error("The size of dest and indices must be the same."); - } - - // src dtype - py::dtype src_dt = src.dtype(); - - // dest dtype - py::dtype dst_dt = dest.dtype(); - - // indices dtype - py::dtype indices_dt = indices.dtype(); - - // validate dtypes - if (dst_dt != src_dt || std::string(indices_dt.str()) != "int32") - { - throw std::runtime_error("bad array dtype"); - } - - - int* indices_data = static_cast(indices_info.ptr); - - std::string src_dt_str = std::string(src_dt.str()); - if (src_dt_str == "float64") - { - double* dest_data = static_cast(dest_info.ptr); - - py::array_t double_array(src); - auto src_data = double_array.unchecked<1>(); - - for (py::ssize_t i = 0; i < dest_info.shape[0]; ++i) - { - dest_data[i] = src_data(static_cast(indices_data[i])); - } - } - else if (src_dt_str == "int32") - { - int* dest_data = static_cast(dest_info.ptr); - - py::array_t int_array(src); - auto src_data = int_array.unchecked<1>(); - - for (py::ssize_t i = 0; i < dest_info.shape[0]; ++i) - { - dest_data[i] = src_data(static_cast(indices_data[i])); - } - } - else if (src_dt_str.rfind(" str_array(src); - //auto src_data = int_array.unchecked<1>(); - - //for (py::ssize_t i = 0; i < dest_info.shape[0]; ++i) - //{ - // dest_data[i] = src_data(static_cast(indices_data[i])); - //} - assert(false); - throw std::runtime_error("str are unsupported in get_value_at_indices"); - } - else - { - assert(false); - throw std::runtime_error("bad array dtype"); - } - - - return dest; -} - -py::array BMIPhreeqcRM::get_value_ptr(std::string name) -{ - if (!this->_initialized) throw NotIntialized(); - - assert(this->var_man); - RMVARS v_enum = this->var_man->GetEnum(name); - if (v_enum == RMVARS::NotFound) throw std::runtime_error("Unknown variable name"); - - BMIVariant& bv = this->var_man->VariantMap[v_enum]; - - if (bv.HasPyArray()) - { - return bv.GetPyArray(); - } - - // this call is REQUIRED for the BMIVariant::GetVoidPtr() to be valid - // and BMIVariant::GetHasPtr() - void* ptr = this->GetValuePtr(name); - - if (!bv.GetHasPtr() && !(bv.GetPType().rfind("language == "Py"); - - //// this call is REQUIRED for the BMIVariant::GetVoidPtr() to be valid - //void* ptr = this->GetValuePtr(name); - - assert(this->GetVarType(name) == bv.GetPType()); - - bool read_only = !(bv.GetHasSetter()); - - py::array a = py::none(); - - if (bv.GetPType() == "float64") - { - assert(bv.GetItemsize() == sizeof(double)); - assert(bv.GetVoidPtr() == this->GetValuePtr(name)); - - py::array_t arr({ bv.GetDim() }, (double*)bv.GetVoidPtr(), py::cast(bv.GetVoidPtr())); - a = arr; - - // Check that the array shares its data - assert(!arr.owndata()); - assert(arr.data() == bv.GetVoidPtr()); - - py::dtype dt = arr.dtype(); - assert(std::string(dt.str()) == bv.GetPType()); - } - else if (bv.GetPType() == "int32") - { - assert(bv.GetItemsize() == sizeof(int)); - assert(bv.GetVoidPtr() == this->GetValuePtr(name)); - - py::array_t arr({ bv.GetDim() }, (int*)bv.GetVoidPtr(), py::cast(bv.GetVoidPtr())); - a = arr; - - // Check that the array shares its data - assert(!arr.owndata()); - assert(arr.data() == bv.GetVoidPtr()); - - py::dtype dt = arr.dtype(); - assert(std::string(dt.str()) == bv.GetPType()); - } - else if (bv.GetPType().rfind(" &strings = bv.GetStringVectorRef(); - assert(strings.size() > 0); - - a = py::cast(strings); - // string arrays are always readonly - read_only = true; - - py::dtype dt = a.dtype(); - assert(std::string(dt.str()) == bv.GetPType()); - } - else - { - assert(false); - } - - if (read_only) - { - assert(a.writeable()); - py::detail::array_proxy(a.ptr())->flags &= ~py::detail::npy_api::NPY_ARRAY_WRITEABLE_; - assert(!a.writeable()); - } - return bv.SetPyArray(a); -} - -void BMIPhreeqcRM::set_value(std::string name, py::array arr) -{ - if (!this->_initialized) throw NotIntialized(); - - RMVARS v_enum = this->var_man->GetEnum(name); - if (v_enum == RMVARS::NotFound) throw std::runtime_error("Unknown variable name"); - - BMIVariant& bv = this->var_man->VariantMap[v_enum]; - - if (!bv.GetInitialized()) - { - this->var_man->task = VarManager::VAR_TASKS::Info; - ((*this->var_man).*bv.GetFn())(); - } - - py::ssize_t ndim = arr.ndim(); - py::dtype dt = arr.dtype(); - std::string dt_str = std::string(dt.str()); - - std::string bv_str = bv.GetPType(); - int dim = bv.GetDim(); - - const py::ssize_t* shape = arr.shape(); - py::ssize_t sz = arr.size(); - - // Check dimension - if (!(ndim > 0 && shape[0] == dim)) - { - std::ostringstream oss; - oss << "dimension error in set_value: " << name; - this->ErrorMessage(oss.str()); - throw std::runtime_error(oss.str()); - } - - if (dt_str == "float64") - { - //std::vector& ref = this->var_man->VarExchange.GetDoubleVectorRef(); - //ref.resize(arr.size()); - //std::memcpy(ref.data(), arr.data(), arr.size() * sizeof(double)); - assert(bv.GetCType() == "double" && dim > 1); - this->SetValue(name, (void*)arr.data()); - } - else if (dt_str == "int32") - { - //std::vector& ref = this->var_man->VarExchange.GetIntVectorRef(); - //ref.resize(arr.size()); - //std::memcpy(ref.data(), arr.data(), arr.size() * sizeof(int)); - assert(bv.GetCType() == "int" && dim > 1); - this->SetValue(name, (void*)arr.data()); - } - else if (dt_str.rfind(" result; - for (auto item : arr) - { - result.push_back(py::cast(item)); - } - - size_t n = result.size(); - assert(n > 0); - assert(shape[0] == n); - assert(ndim == 1); - assert(sz == n); - - if (n == 1) - { - this->SetValue(name, result[0]); - - std::string value; - this->GetValue(name, value); - assert(value == "prefix"); - } - else - { - // this should have been caught in check dims above - assert(false); - } - } - else if (dt_str == "bool") - { - std::vector result; - for (auto item : arr) - { - bool b = py::cast(item); - result.push_back(py::cast(item)); - } - - size_t n = result.size(); - assert(n > 0); - assert(shape[0] == n); - assert(ndim == 1); - assert(sz == n); - - if (n == 1) - { - this->SetValue(name, result[0]); - } - else - { - // this should have been caught in check dims above - assert(false); - } - } - return; -} - -void BMIPhreeqcRM::set_value_at_indices(std::string name, py::array indices, py::array src) -{ - // name and initilized validated in get_value_ptr - py::array dest = get_value_ptr(name); - //{{ - py::buffer_info dest_info = dest.request(); - //}} - - // validate src and indices - py::buffer_info src_info = src.request(); - py::buffer_info indices_info = indices.request(); - - if (src_info.ndim != 1 || indices_info.ndim != 1) - { - throw std::runtime_error("Both src and indices must be 1-dimensional arrays."); - } - - if (src_info.shape[0] != indices_info.shape[0]) - { - throw std::runtime_error("The size of src and indices must be the same."); - } - - // src dtype - py::dtype src_dt = src.dtype(); - - // dest dtype - py::dtype dst_dt = dest.dtype(); - - // indices dtype - py::dtype indices_dt = indices.dtype(); - - // validate dtypes - if (dst_dt != src_dt || std::string(indices_dt.str()) != "int32") - { - throw std::runtime_error("bad array dtype"); - } - - int* indices_data = static_cast(indices_info.ptr); - - std::string src_dt_str = std::string(src_dt.str()); - if (src_dt_str == "float64") - { - double* dest_data = static_cast(dest_info.ptr); - - py::array_t double_array(src); - auto src_data = double_array.unchecked<1>(); - - for (py::ssize_t i = 0; i < src_info.shape[0]; ++i) - { - dest_data[indices_data[i]] = src_data(static_cast(i)); - } - } - else if (src_dt_str == "int32") - { - int* dest_data = static_cast(dest_info.ptr); - - py::array_t int_array(src); - auto src_data = int_array.unchecked<1>(); - - for (py::ssize_t i = 0; i < src_info.shape[0]; ++i) - { - dest_data[indices_data[i]] = src_data(static_cast(i)); - } - } - else if (src_dt_str.rfind(" str_array(src); - //auto src_data = int_array.unchecked<1>(); - - //for (py::ssize_t i = 0; i < dest_info.shape[0]; ++i) - //{ - // dest_data[i] = src_data(static_cast(indices_data[i])); - //} - assert(false); - throw std::runtime_error("str are unsupported in get_value_at_indices"); - } - else - { - assert(false); - throw std::runtime_error("bad array dtype"); - } -} - -PYBIND11_MODULE(phreeqcrm, m) { - - py::class_(m, "bmi_phreeqcrm" , py::dynamic_attr()) - - // Constructor - .def(py::init()) - - // - // Model control functions. - // - - // virtual void Initialize(std::string config_file) = 0; - // def initialize(self, config_file: str) -> None: - .def("initialize", - &BMIPhreeqcRM::Initialize, - initialize_docstring.c_str(), - py::arg("config_file") = py::str("") - ) - - // virtual void Update() = 0; - // def update(self) -> None: - .def("update", - [](BMIPhreeqcRM& self) { - if (!self._initialized) throw NotIntialized(); - self.Update(); - }, - update_docstring.c_str() - ) - - // virtual void UpdateUntil(double time) = 0; - // def update_until(self, time: float) -> None: - .def("update_until", - [](BMIPhreeqcRM& self, double time) { - if (!self._initialized) throw NotIntialized(); - self.UpdateUntil(time); - }, - update_until_docstring.c_str(), - py::arg("time") - ) - - // virtual void Finalize() = 0; - // def finalize(self) -> None: - .def("finalize", - [](BMIPhreeqcRM& self) { - self.Finalize(); - } - ) - - // - // Model information functions. - // - - // virtual std::string GetComponentName() = 0; - // def get_component_name(self) -> str: - .def("get_component_name", - [](BMIPhreeqcRM& self) { - return self.GetComponentName(); - }, - get_component_name_docstring.c_str() - ) - - // virtual int GetInputItemCount() = 0; - // def get_input_item_count(self) -> int: - .def("get_input_item_count", - &BMIPhreeqcRM::GetInputItemCount - ) - - // virtual int GetOutputItemCount() = 0; - // def get_input_item_count(self) -> int: - .def("get_output_item_count", - &BMIPhreeqcRM::GetOutputItemCount - ) - - // virtual std::vector GetInputVarNames() = 0; - // def get_input_var_names(self) -> Tuple[str]: - .def("get_input_var_names", - [](BMIPhreeqcRM& self) { - if (!self._initialized) throw NotIntialized(); - auto output = self.GetInputVarNames(); - py::tuple out = py::cast(output); - return out; - } - ) - - // virtual std::vector GetOutputVarNames() = 0; - // def get_output_var_names(self) -> Tuple[str]: - .def("get_output_var_names", - [](BMIPhreeqcRM& self) { - if (!self._initialized) throw NotIntialized(); - auto output = self.GetOutputVarNames(); - py::tuple out = py::cast(output); - return out; - } - ) - - // - // Variable information functions - // - - // virtual int GetVarGrid(std::string name) = 0; - // def get_var_grid(self, name: str) -> int: - .def("get_var_grid", - &BMIPhreeqcRM::GetVarGrid, - py::arg("name") - ) - - // virtual std::string GetVarType(std::string name) = 0; - // def get_var_type(self, name: str) -> str: - .def("get_var_type", - &BMIPhreeqcRM::GetVarType, - py::arg("name") - ) - - // virtual std::string GetVarUnits(std::string name) = 0; - // def get_var_units(self, name: str) -> str: - .def("get_var_units", - &BMIPhreeqcRM::GetVarUnits, - py::arg("name") - ) - - // virtual int GetVarItemsize(std::string name) = 0; - // def get_var_itemsize(self, name: str) -> int: - .def("get_var_itemsize", - &BMIPhreeqcRM::GetVarItemsize, - py::arg("name") - ) - - - // virtual int GetVarNbytes(std::string name) = 0; - // def get_var_nbytes(self, name: str) -> int: - .def("get_var_nbytes", - &BMIPhreeqcRM::GetVarNbytes, - py::arg("name") - ) - - // virtual std::string GetVarLocation(std::string name) = 0; - // def get_var_location(self, name: str) -> str: - .def("get_var_location", - [](BMIPhreeqcRM& self, std::string name) { - if (!self._initialized) throw NotIntialized(); - return self.GetVarLocation(name); - }, - py::arg("name") - ) - - // virtual double GetCurrentTime() = 0; - // def get_current_time(self) -> float: - .def("get_current_time", - [](BMIPhreeqcRM& self) { - if (!self._initialized) throw NotIntialized(); - return self.GetCurrentTime(); - } - ) - - // virtual double GetStartTime() = 0; - // def get_start_time(self) -> float: - .def("get_start_time", - [](BMIPhreeqcRM& self) { - if (!self._initialized) throw NotIntialized(); - return self.GetStartTime(); - } - ) - - // virtual double GetEndTime() = 0; - // def get_end_time(self) -> float: - .def("get_end_time", - [](BMIPhreeqcRM& self) { - if (!self._initialized) throw NotIntialized(); - return self.GetEndTime(); - } - ) - - // virtual std::string GetTimeUnits() = 0; - // def get_time_units(self) -> str: - .def("get_time_units", - [](BMIPhreeqcRM& self) { - return self.GetTimeUnits(); - } - ) - - // virtual double GetTimeStep() = 0; - // def get_time_step(self) -> float: - .def("get_time_step", - [](BMIPhreeqcRM& self) { - return self.GetTimeStep(); - } - ) - - // - // Variable getters - // - - // virtual void GetValue(std::string name, void *dest) = 0; - // def get_value(self, name: str, dest: np.ndarray) -> np.ndarray: - .def("get_value", - &BMIPhreeqcRM::get_value, - py::arg("name"), - py::arg("dest") - ) - - // virtual void *GetValuePtr(std::string name) = 0; - // def get_value_ptr(self, name: str) -> np.ndarray: - .def("get_value_ptr", - &BMIPhreeqcRM::get_value_ptr, - py::arg("name") - ) - - // virtual void GetValueAtIndices(std::string name, void *dest, int *inds, int count) = 0; - // def get_value_at_indices( - // self, name: str, dest : np.ndarray, inds : np.ndarray - // )->np.ndarray: - .def("get_value_at_indices", - &BMIPhreeqcRM::get_value_at_indices, - py::arg("name"), - py::arg("dest"), - py::arg("inds") - ) - - // - // Variable setters - // - - // virtual void SetValue(std::string name, void *src) = 0; - // def set_value(self, name: str, src: np.ndarray) -> None: - .def("set_value", - &BMIPhreeqcRM::set_value, - py::arg("name"), - py::arg("src") - ) - - // virtual void SetValueAtIndices(std::string name, int *inds, int count, void *src) = 0; - // def set_value_at_indices( - // self, name: str, inds : np.ndarray, src : np.ndarray - // )->None: - .def("set_value_at_indices", - &BMIPhreeqcRM::set_value_at_indices, - py::arg("name"), - py::arg("inds"), - py::arg("src") - ) - - // - // Grid information functions - // - - // virtual int GetGridRank(const int grid) = 0; - // def get_grid_rank(self, grid: int) -> int: - .def("get_grid_rank", - &BMIPhreeqcRM::GetGridRank, - py::arg("grid") - ) - - // virtual int GetGridSize(const int grid) = 0; - // def get_grid_size(self, grid: int) -> int: - .def("get_grid_size", - [](BMIPhreeqcRM& self, int grid) { - if (!self._initialized) throw NotIntialized(); - return self.GetGridSize(grid); - }, - get_grid_size_docstring.c_str(), - py::arg("grid") - ) - - // virtual std::string GetGridType(const int grid) = 0; - // def get_grid_type(self, grid: int) -> str: - .def("get_grid_type", - &BMIPhreeqcRM::GetGridType, - py::arg("grid") - ) - - // void GetGridShape(const int grid, int* shape) override - // def get_grid_shape(self, grid: int, shape: np.ndarray) -> np.ndarray: - .def("get_grid_shape", - [](BMIPhreeqcRM& self, int grid, py::array_t shape) { - throw NotImplemented(); - if (!self._initialized) throw NotIntialized(); - - int* data = shape.mutable_data(); - py::array::ShapeContainer shapeContainer{1}; - shape.reshape(shapeContainer); - data[0] = self.GetGridCellCount(); - return shape; - } - , - //get_grid_shape_docstring.c_str(), - py::arg("grid"), - py::arg("shape") - ) - - // virtual void GetGridSpacing(const int grid, double *spacing) = 0; - // def get_grid_spacing(self, grid: int, spacing: np.ndarray) -> np.ndarray: - .def("get_grid_spacing", - [](BMIPhreeqcRM& self, int grid, py::array_t spacing) { - throw NotImplemented(); - if (!self._initialized) throw NotIntialized(); - return spacing; - } - , - //get_grid_spacing_docstring.c_str(), - py::arg("grid"), - py::arg("spacing") - ) - - // virtual void GetGridOrigin(const int grid, double *origin) = 0; - // def get_grid_origin(self, grid: int, origin: np.ndarray) -> np.ndarray: - .def("get_grid_origin", - [](BMIPhreeqcRM& self, int grid, py::array_t origin) { - throw NotImplemented(); - if (!self._initialized) throw NotIntialized(); - return origin; - } - , - //get_grid_origin_docstring.c_str(), - py::arg("grid"), - py::arg("origin") - ) - - // virtual void GetGridX(const int grid, double* x) = 0; - // def get_grid_x(self, grid: int, x: np.ndarray) -> np.ndarray: - .def("get_grid_x", - [](BMIPhreeqcRM& self, int grid, py::array x) { - throw NotImplemented(); - return x; - }, - //get_grid_x_docstring.c_str(), - py::arg("grid"), - py::arg("x") - ) - - // virtual void GetGridY(const int grid, double* y) = 0; - // def get_grid_y(self, grid: int, y: np.ndarray) -> np.ndarray: - .def("get_grid_y", - [](BMIPhreeqcRM& self, int grid, py::array y) { - throw NotImplemented(); - return y; - }, - //get_grid_y_docstring.c_str(), - py::arg("grid"), - py::arg("y") - ) - - // virtual void GetGridZ(const int grid, double* z) = 0; - // def get_grid_y(self, grid: int, z: np.ndarray) -> np.ndarray: - .def("get_grid_z", - [](BMIPhreeqcRM& self, int grid, py::array z) { - throw NotImplemented(); - return z; - }, - //get_grid_z_docstring.c_str(), - py::arg("grid"), - py::arg("z") - ) - - // virtual int GetGridNodeCount(const int grid) = 0; - // def get_grid_node_count(self, grid: int) -> int: - .def("get_grid_node_count", - [](BMIPhreeqcRM& self, int grid) { - //throw NotImplemented(); - return self.GetGridNodeCount(grid); - }, - // get_grid_node_count_docstring.c_str(), - py::arg("grid") - ) - - // virtual int GetGridEdgeCount(const int grid) = 0; - // def get_grid_edge_count(self, grid: int) -> int: - .def("get_grid_edge_count", - [](BMIPhreeqcRM& self, int grid) { - // throw NotImplemented(); - return self.GetGridEdgeCount(grid); - }, - // get_grid_edge_count_docstring.c_str(), - py::arg("grid") - ) - - // virtual int GetGridFaceCount(const int grid) = 0; - // def get_grid_face_count(self, grid: int) -> int: - .def("get_grid_face_count", - [](BMIPhreeqcRM& self, int grid) { - // throw NotImplemented(); - return self.GetGridFaceCount(grid); - }, - // get_grid_face_count_docstring.c_str(), - py::arg("grid") - ) - - // virtual void GetGridEdgeNodes(const int grid, int* edge_nodes) = 0; - // def get_grid_edge_nodes(self, grid: int, edge_nodes: np.ndarray) -> np.ndarray: - .def("get_grid_edge_nodes", - [](BMIPhreeqcRM& self, int grid, py::array_t edge_nodes) { - throw NotImplemented(); - return edge_nodes; - }, - // get_grid_edge_nodes_docstring.c_str(), - py::arg("grid"), - py::arg("edge_nodes") - ) - - // virtual void GetGridFaceEdges(const int grid, int* face_edges) = 0; - // def get_grid_face_edges(self, grid: int, face_edges: np.ndarray) -> np.ndarray: - .def("get_grid_face_edges", - [](BMIPhreeqcRM& self, int grid, py::array_t face_edges) { - throw NotImplemented(); - return face_edges; - }, - // get_grid_face_edges_docstring.c_str(), - py::arg("grid"), - py::arg("face_edges") - ) - - // virtual void GetGridFaceNodes(const int grid, int* face_nodes) = 0; - // def get_grid_face_nodes(self, grid: int, face_nodes: np.ndarray) -> np.ndarray: - .def("get_grid_face_nodes", - [](BMIPhreeqcRM& self, int grid, py::array_t face_nodes) { - throw NotImplemented(); - return face_nodes; - }, - // get_grid_face_nodes_docstring.c_str(), - py::arg("grid"), - py::arg("face_nodes") - ) - - // virtual void GetGridNodesPerFace(const int grid, int* nodes_per_face) = 0; - // def get_grid_nodes_per_face( - // self, grid: int, nodes_per_face: np.ndarray - // ) -> np.ndarray: - .def("get_grid_nodes_per_face", - [](BMIPhreeqcRM& self, int grid, py::array_t nodes_per_face) { - throw NotImplemented(); - return nodes_per_face; - }, - // get_grid_nodes_per_face_docstring.c_str(), - py::arg("grid"), - py::arg("nodes_per_face") - ) - - // - // Extras - // - - - .def("get_components", - [](BMIPhreeqcRM& self) { - // initialization not necessary - auto output = self.GetComponents(); - py::tuple out = py::cast(output); - return out; - }, - get_components_docstring.c_str() - ) - - // int GetGridCellCount(void) - // def GetGridCellCount(self: phreeqcrm.bmi_phreeqcrm) -> int: - .def("GetGridCellCount", - &BMIPhreeqcRM::GetGridCellCount - ) - - // int GetMpiMyself(void) const - // def GetMpiMyself(self: phreeqcrm.bmi_phreeqcrm) -> int: - .def("GetMpiMyself", - &BMIPhreeqcRM::GetMpiMyself - ) - - // int GetThreadCount(void) - // def GetThreadCount(self: phreeqcrm.bmi_phreeqcrm) -> int: - .def("GetThreadCount", - &BMIPhreeqcRM::GetThreadCount - ) - - - // IRM_RESULT SetComponentH2O(bool tf); - // SetComponentH2O(self: phreeqcrm.bmi_phreeqcrm, arg0: bool) -> IRM_RESULT - .def("SetComponentH2O", - &BMIPhreeqcRM::SetComponentH2O, - py::arg("tf") - ) - - .def_readonly("_initialized", &BMIPhreeqcRM::_initialized) - ; - - py::enum_(m, "IRM_RESULT") - .value("IRM_OK", IRM_OK) - .value("IRM_OUTOFMEMORY", IRM_OUTOFMEMORY) - .value("IRM_BADVARTYPE", IRM_BADVARTYPE) - .value("IRM_INVALIDARG", IRM_INVALIDARG) - .value("IRM_INVALIDROW", IRM_INVALIDROW) - .value("IRM_INVALIDCOL", IRM_INVALIDCOL) - .value("IRM_BADINSTANCE", IRM_BADINSTANCE) - .value("IRM_FAIL", IRM_FAIL) - .export_values(); -} diff --git a/pybind/pybind11 b/pybind/pybind11 deleted file mode 160000 index be97c5a98..000000000 --- a/pybind/pybind11 +++ /dev/null @@ -1 +0,0 @@ -Subproject commit be97c5a98b4b252c524566f508b5c79410d118c6 diff --git a/pybind/tests/AdvectBMI_py.yaml b/pybind/tests/AdvectBMI_py.yaml deleted file mode 100644 index 65cc2085c..000000000 --- a/pybind/tests/AdvectBMI_py.yaml +++ /dev/null @@ -1,1163 +0,0 @@ -- key: SetGridCellCount - count: 40 -- key: ThreadCount - nthreads: 3 -- key: SetErrorHandlerMode - mode: 1 -- key: SetComponentH2O - tf: false -- key: SetRebalanceFraction - f: 0.5 -- key: SetRebalanceByCell - tf: true -- key: UseSolutionDensityVolume - tf: false -- key: SetPartitionUZSolids - tf: false -- key: SetFilePrefix - prefix: AdvectBMI_py -- key: OpenFiles -- key: SetUnitsSolution - option: 2 -- key: SetUnitsPPassemblage - option: 1 -- key: SetUnitsExchange - option: 1 -- key: SetUnitsSurface - option: 1 -- key: SetUnitsGasPhase - option: 1 -- key: SetUnitsSSassemblage - option: 1 -- key: SetUnitsKinetics - option: 1 -- key: SetTimeConversion - conv_factor: 1.1574074074074073e-05 -- key: SetRepresentativeVolume - rv: - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 -- key: SetDensityUser - density: - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 - - 1.0 -- key: SetPorosity - por: - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 - - 0.2 -- key: SetSaturationUser - sat: - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 -- key: SetPrintChemistryMask - cell_mask: - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 0 - - 0 - - 0 - - 0 - - 0 - - 0 - - 0 - - 0 - - 0 - - 0 - - 0 - - 0 - - 0 - - 0 - - 0 - - 0 - - 0 - - 0 - - 0 - - 0 -- key: CreateMapping - grid2chem: - - 0 - - 1 - - 2 - - 3 - - 4 - - 5 - - 6 - - 7 - - 8 - - 9 - - 10 - - 11 - - 12 - - 13 - - 14 - - 15 - - 16 - - 17 - - 18 - - 19 - - 0 - - 1 - - 2 - - 3 - - 4 - - 5 - - 6 - - 7 - - 8 - - 9 - - 10 - - 11 - - 12 - - 13 - - 14 - - 15 - - 16 - - 17 - - 18 - - 19 -- key: SetPrintChemistryOn - workers: false - initial_phreeqc: true - utility: false -- key: LoadDatabase - database: phreeqc.dat -- key: RunFile - workers: true - initial_phreeqc: true - utility: true - chemistry_name: advect.pqi -- key: RunString - workers: true - initial_phreeqc: false - utility: true - input_string: DELETE; -all -- key: AddOutputVars - option: AddOutputVars - definition: 'true' -- key: FindComponents -- key: InitialPhreeqc2Module_mix - ic1: - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - ic2: - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - - -1 - f1: - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 - - 1 -- key: InitialPhreeqcCell2Module - n: -1 - cell_numbers: - - 18 - - 19 -- key: SetTimeStep - time_step: 0.0 -- key: SetTime - time: 0.0 -- key: RunCells -- key: SetTimeStep - time_step: 86400.0 diff --git a/pybind/tests/advect.pqi b/pybind/tests/advect.pqi deleted file mode 100644 index 56e711bf4..000000000 --- a/pybind/tests/advect.pqi +++ /dev/null @@ -1,31 +0,0 @@ -TITLE Example 11.--Transport and cation exchange. -SOLUTION 0 CaCl2 - units mmol/kgw - temp 25.0 - pH 7.0 charge - pe 12.5 O2(g) -0.68 - Ca 0.6 - Cl 1.2 -SOLUTION 1-40 Initial solution for column - units mmol/kgw - temp 25.0 - pH 7.0 charge - pe 12.5 O2(g) -0.68 - Na 1.0 - K 0.2 - N(5) 1.2 -END -EXCHANGE 1-40 - -equilibrate 1 - X 0.0011 -END -SELECTED_OUTPUT - -file advect.sel - -reset false -# -step - -totals Na Cl K Ca -USER_PUNCH - -heading Temperature Pressure Hyd_K - 10 PUNCH TC, PRESSURE - 20 PUNCH CALLBACK(cell_no, 0, "HYDRAULIC_K") -END diff --git a/pybind/tests/constants.py b/pybind/tests/constants.py deleted file mode 100644 index 3dd3f0f3e..000000000 --- a/pybind/tests/constants.py +++ /dev/null @@ -1,19 +0,0 @@ -# constants.py -class Constants: - def __init__(self): - file_directory = os.path.dirname(__file__) - self._yaml = os.path.join(file_directory, 'AdvectBMI_py.yaml') - self._database = os.path.join(file_directory, 'phreeqc.dat') - self._pqi = os.path.join(file_directory, 'advect_pqi') - - @property - def yaml(self): - return self._yaml - - @property - def database(self): - return self._database - - @property - def pqi(self): - return self._pqi diff --git a/pybind/tests/phreeqc.dat b/pybind/tests/phreeqc.dat deleted file mode 100644 index 4883650d0..000000000 --- a/pybind/tests/phreeqc.dat +++ /dev/null @@ -1,1930 +0,0 @@ -# PHREEQC.DAT for calculating temperature and pressure dependence of reactions, and the specific conductance and viscosity of the solution. Based on: -# diffusion coefficients and molal volumina of aqueous species, solubility and volume of minerals, and critical temperatures and pressures of gases in Peng-Robinson's EOS. -# Details are given at the end of this file. - -SOLUTION_MASTER_SPECIES -# -#element species alk gfw_formula element_gfw -# -H H+ -1.0 H 1.008 -H(0) H2 0 H -H(1) H+ -1.0 0 -E e- 0 0.0 0 -O H2O 0 O 16.0 -O(0) O2 0 O -O(-2) H2O 0 0 -Ca Ca+2 0 Ca 40.08 -Mg Mg+2 0 Mg 24.312 -Na Na+ 0 Na 22.9898 -K K+ 0 K 39.102 -Fe Fe+2 0 Fe 55.847 -Fe(+2) Fe+2 0 Fe -Fe(+3) Fe+3 -2.0 Fe -Mn Mn+2 0 Mn 54.938 -Mn(+2) Mn+2 0 Mn -Mn(+3) Mn+3 0 Mn -Al Al+3 0 Al 26.9815 -Ba Ba+2 0 Ba 137.34 -Sr Sr+2 0 Sr 87.62 -Si H4SiO4 0 SiO2 28.0843 -Cl Cl- 0 Cl 35.453 -C CO3-2 2.0 HCO3 12.0111 -C(+4) CO3-2 2.0 HCO3 -C(-4) CH4 0 CH4 -Alkalinity CO3-2 1.0 Ca0.5(CO3)0.5 50.05 -S SO4-2 0 SO4 32.064 -S(6) SO4-2 0 SO4 -S(-2) HS- 1.0 S -N NO3- 0 N 14.0067 -N(+5) NO3- 0 N -N(+3) NO2- 0 N -N(0) N2 0 N -N(-3) NH4+ 0 N 14.0067 -#Amm AmmH+ 0 AmmH 17.031 -B H3BO3 0 B 10.81 -P PO4-3 2.0 P 30.9738 -F F- 0 F 18.9984 -Li Li+ 0 Li 6.939 -Br Br- 0 Br 79.904 -Zn Zn+2 0 Zn 65.37 -Cd Cd+2 0 Cd 112.4 -Pb Pb+2 0 Pb 207.19 -Cu Cu+2 0 Cu 63.546 -Cu(+2) Cu+2 0 Cu -Cu(+1) Cu+1 0 Cu -# redox-uncoupled gases -Hdg Hdg 0 Hdg 2.016 # H2 gas -Oxg Oxg 0 Oxg 32 # O2 gas -Mtg Mtg 0 Mtg 16.032 # CH4 gas -Sg H2Sg 1.0 H2Sg 34.08 -Ntg Ntg 0 Ntg 28.0134 # N2 gas - -SOLUTION_SPECIES -H+ = H+ - -gamma 9.0 0 - -dw 9.31e-9 1000 0.46 1e-10 # The dw parameters are defined in ref. 3. -# Dw(TK) = 9.31e-9 * exp(1000 / TK - 1000 / 298.15) * viscos_0_25 / viscos_0_tc -# Dw(I) = Dw(TK) * exp(-0.46 * DH_A * |z_H+| * I^0.5 / (1 + DH_B * I^0.5 * 1e-10 / (1 + I^0.75))) - -viscosity 9.35e-2 -8.31e-2 2.487e-2 4.49e-4 2.01e-2 1.570 # for viscosity parameters see ref. 4 -e- = e- -H2O = H2O -# H2O + 0.01e- = H2O-0.01; -log_k -9 # aids convergence -Ca+2 = Ca+2 - -gamma 5.0 0.1650 - -dw 0.793e-9 97 3.4 24.6 - -Vm -0.3456 -7.252 6.149 -2.479 1.239 5 1.60 -57.1 -6.12e-3 1 # The apparent volume parameters are defined in ref. 1 & 2 - -viscosity 0.359 -0.158 4.2e-2 1.5e-3 8.04e-3 2.30 # ref. 4, CaCl2 < 6 M -Mg+2 = Mg+2 - -gamma 5.5 0.20 - -dw 0.705e-9 111 2.4 13.7 - -Vm -1.410 -8.6 11.13 -2.39 1.332 5.5 1.29 -32.9 -5.86e-3 1 - -viscosity 0.426 0 0 1.66e-3 4.32e-3 2.461 -Na+ = Na+ - -gamma 4.0 0.075 - -gamma 4.08 0.082 # halite solubility - -dw 1.33e-9 122 1.52 3.70 - -Vm 2.28 -4.38 -4.1 -0.586 0.09 4 0.3 52 -3.33e-3 0.566 -# for calculating densities (rho) when I > 3... - # -Vm 2.28 -4.38 -4.1 -0.586 0.09 4 0.3 52 -3.33e-3 0.45 - -viscosity 0.1387 -8.66e-2 1.25e-2 1.45e-2 7.5e-3 1.062 -K+ = K+ - -gamma 3.5 0.015 - -dw 1.96e-9 395 2.5 21 - -Vm 3.322 -1.473 6.534 -2.712 9.06e-2 3.5 0 29.7 0 1 - -viscosity 0.116 -0.191 1.52e-2 1.40e-2 2.59e-2 0.9028 -Fe+2 = Fe+2 - -gamma 6.0 0 - -dw 0.719e-9 - -Vm -0.3255 -9.687 1.536 -2.379 0.3033 6 -4.21e-2 39.7 0 1 -Mn+2 = Mn+2 - -gamma 6.0 0 - -dw 0.688e-9 - -Vm -1.10 -8.03 4.08 -2.45 1.4 6 8.07 0 -1.51e-2 0.118 -Al+3 = Al+3 - -gamma 9.0 0 - -dw 0.559e-9 - -Vm -2.28 -17.1 10.9 -2.07 2.87 9 0 0 5.5e-3 1 # ref. 2 and Barta and Hepler, 1986, Can. J.C. 64, 353. -Ba+2 = Ba+2 - -gamma 5.0 0 - -gamma 4.0 0.153 # Barite solubility - -dw 0.848e-9 100 - -Vm 2.063 -10.06 1.9534 -2.36 0.4218 5 1.58 -12.03 -8.35e-3 1 - -viscosity 0.338 -0.227 1.39e-2 3.07e-2 0 0.768 -Sr+2 = Sr+2 - -gamma 5.260 0.121 - -dw 0.794e-9 161 - -Vm -1.57e-2 -10.15 10.18 -2.36 0.860 5.26 0.859 -27.0 -4.1e-3 1.97 - -viscosity 0.472 -0.252 5.51e-3 3.67e-3 0 1.876 -H4SiO4 = H4SiO4 - -dw 1.10e-9 - -Vm 10.5 1.7 20 -2.7 0.1291 # supcrt + 2*H2O in a1 -Cl- = Cl- - -gamma 3.5 0.015 - -gamma 3.63 0.017 # cf. pitzer.dat - -dw 2.03e-9 194 1.6 6.9 - -Vm 4.465 4.801 4.325 -2.847 1.748 0 -0.331 20.16 0 1 - -viscosity 0 0 0 0 0 0 1 # the reference solute -CO3-2 = CO3-2 - -gamma 5.4 0 - -dw 0.955e-9 27.4 13.7 94.1 - -Vm 8.61 -10.26 -19.54 -0.150 4.63 0 3.32 0 -3.56e-2 0.770 - -viscosity 0 0.289 3.70e-2 5e-5 -3.03e-2 2.013 -2.04 -SO4-2 = SO4-2 - -gamma 5.0 -0.04 - -dw 1.07e-9 187 2.64 22.6 - -Vm 9.379 3.26 0 -7.13 4.30 0 0 0 -3.73e-2 0 # with analytical_expressions for log K of NaSO4-, KSO4- & MgSO4, 0 - 200 oC - -viscosity -1.83 1.907 4.8e-4 1.7e-3 -1.60e-2 4.40 -0.143 -NO3- = NO3- - -gamma 3.0 0 - -dw 1.9e-9 184 1.85 3.85 - -Vm 6.32 6.78 0 -3.06 0.346 0 0.93 0 -0.012 1 - -viscosity 8.37e-2 -0.458 1.54e-2 0.340 1.79e-2 5.02e-2 0.7381 -#AmmH+ = AmmH+ -# -gamma 2.5 0 -# -dw 1.98e-9 312 0.95 4.53 -# -Vm 4.837 2.345 5.522 -2.88 1.096 3 -1.456 75.0 7.17e-3 1 -# -viscosity 7.25e-2 -0.142 1.97e-2 8.44e-3 3.92e-2 0.945 -H3BO3 = H3BO3 - -dw 1.1e-9 - -Vm 7.0643 8.8547 3.5844 -3.1451 -.2000 # supcrt -PO4-3 = PO4-3 - -gamma 4.0 0 - -dw 0.612e-9 - -Vm 1.24 -9.07 9.31 -2.4 5.61 0 0 0 -1.41e-2 1 -F- = F- - -gamma 3.5 0 - -dw 1.46e-9 10 - -Vm 0.928 1.36 6.27 -2.84 1.84 0 0 -0.318 0 1 -Li+ = Li+ - -gamma 6.0 0 - -dw 1.03e-9 80 - -Vm -0.419 -0.069 13.16 -2.78 0.416 0 0.296 -12.4 -2.74e-3 1.26 # ref. 2 and Ellis, 1968, J. Chem. Soc. A, 1138 - -viscosity 0.162 -2.45e-2 3.73e-2 9.7e-4 8.1e-4 2.087 -Br- = Br- - -gamma 3.0 0 - -dw 2.01e-9 258 - -Vm 6.72 2.85 4.21 -3.14 1.38 0 -9.56e-2 7.08 -1.56e-3 1 - -viscosity -1.15e-2 -5.75e-2 5.72e-2 1.46e-2 0.116 0.9295 0.820 -Zn+2 = Zn+2 - -gamma 5.0 0 - -dw 0.715e-9 - -Vm -1.96 -10.4 14.3 -2.35 1.46 5 -1.43 24 1.67e-2 1.11 -Cd+2 = Cd+2 - -dw 0.717e-9 - -Vm 1.63 -10.7 1.01 -2.34 1.47 5 0 0 0 1 -Pb+2 = Pb+2 - -dw 0.945e-9 - -Vm -.0051 -7.7939 8.8134 -2.4568 1.0788 4.5 # supcrt -Cu+2 = Cu+2 - -gamma 6.0 0 - -dw 0.733e-9 - -Vm -1.13 -10.5 7.29 -2.35 1.61 6 9.78e-2 0 3.42e-3 1 -# redox-uncoupled gases -Hdg = Hdg # H2 - -dw 5.13e-9 - -Vm 6.52 0.78 0.12 # supcrt -Oxg = Oxg # O2 - -dw 2.35e-9 - -Vm 5.7889 6.3536 3.2528 -3.0417 -0.3943 # supcrt -Mtg = Mtg # CH4 - -dw 1.85e-9 - -Vm 9.01 -1.11 0 -1.85 -1.50 # Hnedkovsky et al., 1996, JCT 28, 125 -Ntg = Ntg # N2 - -dw 1.96e-9 - -Vm 7 # Pray et al., 1952, IEC 44. 1146 -H2Sg = H2Sg # H2S - -dw 2.1e-9 - -Vm 1.39 28.3 0 -7.22 -0.59 # Hnedkovsky et al., 1996, JCT 28, 125 -# aqueous species -H2O = OH- + H+ - -analytic 293.29227 0.1360833 -10576.913 -123.73158 0 -6.996455e-5 - -gamma 3.5 0 - -dw 5.27e-9 548 0.52 1e-10 - -Vm -9.66 28.5 80.0 -22.9 1.89 0 1.09 0 0 1 - -viscosity -1.02e-1 0.189 9.4e-3 -4e-5 0 3.281 -2.053 # < 5 M Li,Na,KOH -2 H2O = O2 + 4 H+ + 4 e- - -log_k -86.08 - -delta_h 134.79 kcal - -dw 2.35e-9 - -Vm 5.7889 6.3536 3.2528 -3.0417 -0.3943 # supcrt -2 H+ + 2 e- = H2 - -log_k -3.15 - -delta_h -1.759 kcal - -dw 5.13e-9 - -Vm 6.52 0.78 0.12 # supcrt -H+ + Cl- = HCl - -log_k -0.5 - -analytical_expression 0.334 -2.684e-3 1.015 # from Pitzer.dat, up to 15 M HCl, 0 - 50°C - -gamma 0 0.4256 - -viscosity 0.921 -0.765 8.32e-3 8.25e-4 2.53e-3 4.223 -CO3-2 + H+ = HCO3- - -log_k 10.329 - -delta_h -3.561 kcal - -analytic 107.8871 0.03252849 -5151.79 -38.92561 563713.9 - -gamma 5.4 0 - -dw 1.18e-9 -163 0.808 -3.18 - -Vm 9.14 -1.64 -12.00 0 1.63 0 0 132 0 0.667 - -viscosity 0 0.670 1.03e-2 0 0 0 1.082 -CO3-2 + 2 H+ = CO2 + H2O - -log_k 16.681 - -delta_h -5.738 kcal - -analytic 464.1965 0.09344813 -26986.16 -165.75951 2248628.9 - -dw 1.92e-9 - -Vm 7.29 0.92 2.07 -1.23 -1.60 # McBride et al. 2015, JCED 60, 171 -2CO2 = (CO2)2 # activity correction for CO2 solubility at high P, T - -log_k -1.8 - -analytical_expression 8.68 -0.0103 -2190 - -Vm 14.58 1.84 4.14 -2.46 -3.20 -CO3-2 + 10 H+ + 8 e- = CH4 + 3 H2O - -log_k 41.071 - -delta_h -61.039 kcal - -dw 1.85e-9 - -Vm 9.01 -1.11 0 -1.85 -1.50 # Hnedkovsky et al., 1996, JCT 28, 125 -SO4-2 + H+ = HSO4- - -log_k 1.988 - -delta_h 3.85 kcal - -analytic -56.889 0.006473 2307.9 19.8858 - -dw 1.33e-9 - -Vm 8.2 9.2590 2.1108 -3.1618 1.1748 0 -0.3 15 0 1 -HS- = S-2 + H+ - -log_k -12.918 - -delta_h 12.1 kcal - -gamma 5.0 0 - -dw 0.731e-9 -SO4-2 + 9 H+ + 8 e- = HS- + 4 H2O - -log_k 33.65 - -delta_h -60.140 kcal - -gamma 3.5 0 - -dw 1.73e-9 - -Vm 5.0119 4.9799 3.4765 -2.9849 1.4410 # supcrt -HS- + H+ = H2S - -log_k 6.994 - -delta_h -5.30 kcal - -analytical -11.17 0.02386 3279.0 - -dw 2.1e-9 - -Vm 1.39 28.3 0 -7.22 -0.59 # Hnedkovsky et al., 1996, JCT 28, 125 -2H2S = (H2S)2 # activity correction for H2S solubility at high P, T - -analytical_expression 10.227 -0.01384 -2200 - -Vm 36.41 -71.95 0 0 2.58 -H2Sg = HSg- + H+ - -log_k -6.994 - -delta_h 5.30 kcal - -analytical_expression 11.17 -0.02386 -3279.0 - -gamma 3.5 0 - -dw 1.73e-9 - -Vm 5.0119 4.9799 3.4765 -2.9849 1.4410 # supcrt -2H2Sg = (H2Sg)2 # activity correction for H2S solubility at high P, T - -analytical_expression 10.227 -0.01384 -2200 - -Vm 36.41 -71.95 0 0 2.58 -NO3- + 2 H+ + 2 e- = NO2- + H2O - -log_k 28.570 - -delta_h -43.760 kcal - -gamma 3.0 0 - -dw 1.91e-9 - -Vm 5.5864 5.8590 3.4472 -3.0212 1.1847 # supcrt -2 NO3- + 12 H+ + 10 e- = N2 + 6 H2O - -log_k 207.08 - -delta_h -312.130 kcal - -dw 1.96e-9 - -Vm 7 # Pray et al., 1952, IEC 44. 1146 -#AmmH+ = Amm + H+ -NO3- + 10 H+ + 8 e- = NH4+ + 3 H2O - -log_k 119.077 - -delta_h -187.055 kcal - -gamma 2.5 0 - -dw 1.98e-9 312 0.95 4.53 - -Vm 4.837 2.345 5.522 -2.88 1.096 3 -1.456 75.0 7.17e-3 1 - -viscosity 7.25e-2 -0.142 1.97e-2 8.44e-3 3.92e-2 0.945 - -NH4+ = NH3 + H+ - -log_k -9.252 - -delta_h 12.48 kcal - -analytic 0.6322 -0.001225 -2835.76 - -dw 2.28e-9 - -Vm 6.69 2.8 3.58 -2.88 1.43 -#NO3- + 10 H+ + 8 e- = AmmH+ + 3 H2O -# -log_k 119.077 -# -delta_h -187.055 kcal -# -gamma 2.5 0 -# -Vm 4.837 2.345 5.522 -2.88 1.096 3 -1.456 75.0 7.17e-3 1 -#AmmH+ + SO4-2 = AmmHSO4- -NH4+ + SO4-2 = NH4SO4- - -log_k 1.1 - -delta_h -0.47 kcal - -gamma 0 0 - -Vm 13.69 0 -33.54 0 0 0 11.99 0 -0.134 1 - -dw 7.46e-10 - -viscosity -0.109 0.242 1.218e-3 -3.14e-2 8.9e-3 1.631 0.255 -H3BO3 = H2BO3- + H+ - -log_k -9.24 - -delta_h 3.224 kcal -H3BO3 + F- = BF(OH)3- - -log_k -0.4 - -delta_h 1.850 kcal -H3BO3 + 2 F- + H+ = BF2(OH)2- + H2O - -log_k 7.63 - -delta_h 1.618 kcal -H3BO3 + 2 H+ + 3 F- = BF3OH- + 2 H2O - -log_k 13.67 - -delta_h -1.614 kcal -H3BO3 + 3 H+ + 4 F- = BF4- + 3 H2O - -log_k 20.28 - -delta_h -1.846 kcal -PO4-3 + H+ = HPO4-2 - -log_k 12.346 - -delta_h -3.530 kcal - -gamma 5.0 0 - -dw 0.69e-9 - -Vm 3.52 1.09 8.39 -2.82 3.34 0 0 0 0 1 -PO4-3 + 2 H+ = H2PO4- - -log_k 19.553 - -delta_h -4.520 kcal - -gamma 5.4 0 - -dw 0.846e-9 - -Vm 5.58 8.06 12.2 -3.11 1.3 0 0 0 1.62e-2 1 -PO4-3 + 3H+ = H3PO4 - log_k 21.721 # log_k and delta_h from minteq.v4.dat, NIST46.3 - delta_h -10.1 kJ - -Vm 7.47 12.4 6.29 -3.29 0 -H+ + F- = HF - -log_k 3.18 - -delta_h 3.18 kcal - -analytic -2.033 0.012645 429.01 - -Vm 3.4753 .7042 5.4732 -2.8081 -.0007 # supcrt -H+ + 2 F- = HF2- - -log_k 3.76 - -delta_h 4.550 kcal - -Vm 5.2263 4.9797 3.7928 -2.9849 1.2934 # supcrt -Ca+2 + H2O = CaOH+ + H+ - -log_k -12.78 -Ca+2 + CO3-2 = CaCO3 - -log_k 3.224 - -delta_h 3.545 kcal - -analytic -1228.732 -0.299440 35512.75 485.818 - -dw 4.46e-10 # complexes: calc'd with the Pikal formula - -Vm -.2430 -8.3748 9.0417 -2.4328 -.0300 # supcrt -Ca+2 + CO3-2 + H+ = CaHCO3+ - -log_k 11.435 - -delta_h -0.871 kcal - -analytic 1317.0071 0.34546894 -39916.84 -517.70761 563713.9 - -gamma 6.0 0 - -dw 5.06e-10 - -Vm 3.1911 .0104 5.7459 -2.7794 .3084 5.4 # supcrt -Ca+2 + SO4-2 = CaSO4 - -log_k 2.25 - -delta_h 1.325 kcal - -dw 4.71e-10 - -Vm 2.7910 -.9666 6.1300 -2.7390 -.0010 # supcrt -Ca+2 + HSO4- = CaHSO4+ - -log_k 1.08 -Ca+2 + PO4-3 = CaPO4- - -log_k 6.459 - -delta_h 3.10 kcal - -gamma 5.4 0.0 -Ca+2 + HPO4-2 = CaHPO4 - -log_k 2.739 - -delta_h 3.3 kcal -Ca+2 + H2PO4- = CaH2PO4+ - -log_k 1.408 - -delta_h 3.4 kcal - -gamma 5.4 0.0 -# Ca+2 + F- = CaF+ - # -log_k 0.94 - # -delta_h 4.120 kcal - # -gamma 5.5 0.0 - # -Vm .9846 -5.3773 7.8635 -2.5567 .6911 5.5 # supcrt -Mg+2 + H2O = MgOH+ + H+ - -log_k -11.44 - -delta_h 15.952 kcal - -gamma 6.5 0 -Mg+2 + CO3-2 = MgCO3 - -log_k 2.98 - -delta_h 2.713 kcal - -analytic 0.9910 0.00667 - -dw 4.21e-10 - -Vm -.5837 -9.2067 9.3687 -2.3984 -.0300 # supcrt -Mg+2 + H+ + CO3-2 = MgHCO3+ - -log_k 11.399 - -delta_h -2.771 kcal - -analytic 48.6721 0.03252849 -2614.335 -18.00263 563713.9 - -gamma 4.0 0 - -dw 4.78e-10 - -Vm 2.7171 -1.1469 6.2008 -2.7316 .5985 4 # supcrt -Mg+2 + SO4-2 = MgSO4 - -log_k 2.42; -delta_h 19.0 kJ - -analytical_expression 0 9.64e-3 -136 # mean salt gamma from Pitzer.dat and epsomite/hexahydrite/kieserite solubilities, 0 - 200 oC - -gamma 0 0.20 - -Vm 13.18 -25.67 -21.23 0 0.800 0 0 0 0 0 - -dw 4.45e-10 - -viscosity -0.590 0.768 -3.8e-4 0.283 1.1e-3 1.09 0 -SO4-2 + MgSO4 = Mg(SO4)2-2 - -log_k 0.52; -delta_h -13.6 kJ - -analytical_expression 0 -1.51e-3 0 0 8.604e4 # mean salt gamma from Pitzer.dat and epsomite/hexahydrite/kieserite solubilities, 0 - 200 oC - -gamma 7 0.047 - -Vm 12.725 -28.73 0.219 0 -0.264 0 23.44 0 0.213 5.1e-2 - -Dw 1e-9 -2926 6.10e-2 -5.41 - -viscosity -0.162 9.6e-4 -4.65e-2 0.179 1.56e-2 1.66 0 -Mg+2 + PO4-3 = MgPO4- - -log_k 6.589 - -delta_h 3.10 kcal - -gamma 5.4 0 -Mg+2 + HPO4-2 = MgHPO4 - -log_k 2.87 - -delta_h 3.3 kcal -Mg+2 + H2PO4- = MgH2PO4+ - -log_k 1.513 - -delta_h 3.4 kcal - -gamma 5.4 0 -Mg+2 + F- = MgF+ - -log_k 1.82 - -delta_h 3.20 kcal - -gamma 4.5 0 - -Vm .6494 -6.1958 8.1852 -2.5229 .9706 4.5 # supcrt -Na+ + OH- = NaOH - -log_k -10 # remove this complex -# Na+ + CO3-2 = NaCO3- # the HCO3- and CO3-2 cmplxs are not necessary for the SC - # -log_k 1.27 - # -delta_h 8.91 kcal - # -dw 1.2e-9 -400 1e-10 1e-10 - # -Vm 3.812 0.196 20.0 -9.60 3.02 1e-5 2.65 0 2.54e-2 1 - # -viscosity 0.104 -1.65 0.169 8.66e-2 2.60e-2 1.76 -0.90 -# Na+ + HCO3- = NaHCO3 - # -log_k 0.14 - # -delta_h -6.71 kcal - # -dw 6.73e-10 -400 1e-10 1e-10 - # -Vm 6.22 - # -viscosity -0.026 0 0 -0.182 0 3 -Na+ + SO4-2 = NaSO4- - -log_k 0.6; -delta_h -14.4 kJ - -analytical_expression -7.99 1.637e-2 0 0 3.29e5 # mirabilite/thenardite solubilities, 0 - 200 oC - -gamma 0 0 - -Vm 9.993 -8.75 0 -2.95 2.59 0 8.40 0 -1.82e-2 0.672 - -dw 1.183e-9 438 1e-10 1e-10 - -viscosity 7.94e-2 6.96e-2 1.51e-2 7.62e-2 2.84e-2 1.74 0.120 -Na+ + HPO4-2 = NaHPO4- - -log_k 0.29 - -gamma 5.4 0 - -Vm 5.2 8.1 13 -3 0.9 0 0 1.62e-2 1 -Na+ + F- = NaF - -log_k -0.24 - -Vm 2.7483 -1.0708 6.1709 -2.7347 -.030 # supcrt -K+ + SO4-2 = KSO4- - -log_k 0.6; -delta_h -10.4 kJ - -analytical_expression -4.022 8.217e-3 0 0 1.90e5 # arcanite solubility, 0 - 200 oC - -gamma 0 8.3e-3 - -Vm 8.942 -5.05 -15.03 0 3.61 0 25.14 0 -5.06e-2 0.166 - -dw 5.11e-10 1694 -0.587 -4.43 - -viscosity -2.71 3.09 6e-4 -0.629 9.38e-2 0.778 0.975 -K+ + HPO4-2 = KHPO4- - -log_k 0.29 - -gamma 5.4 0 - -Vm 5.4 8.1 19 -3.1 0.7 0 0 0 1.62e-2 1 -Fe+2 + H2O = FeOH+ + H+ - -log_k -9.5 - -delta_h 13.20 kcal - -gamma 5.0 0 -Fe+2 + 3H2O = Fe(OH)3- + 3H+ - -log_k -31.0 - -delta_h 30.3 kcal - -gamma 5.0 0 -Fe+2 + Cl- = FeCl+ - -log_k 0.14 -Fe+2 + CO3-2 = FeCO3 - -log_k 4.38 -Fe+2 + HCO3- = FeHCO3+ - -log_k 2.0 -Fe+2 + SO4-2 = FeSO4 - -log_k 2.25 - -delta_h 3.230 kcal - -Vm -13 0 123 -Fe+2 + HSO4- = FeHSO4+ - -log_k 1.08 -Fe+2 + 2HS- = Fe(HS)2 - -log_k 8.95 -Fe+2 + 3HS- = Fe(HS)3- - -log_k 10.987 -Fe+2 + HPO4-2 = FeHPO4 - -log_k 3.6 -Fe+2 + H2PO4- = FeH2PO4+ - -log_k 2.7 - -gamma 5.4 0 -Fe+2 + F- = FeF+ - -log_k 1.0 -Fe+2 = Fe+3 + e- - -log_k -13.02 - -delta_h 9.680 kcal - -gamma 9.0 0 -Fe+3 + H2O = FeOH+2 + H+ - -log_k -2.19 - -delta_h 10.4 kcal - -gamma 5.0 0 -Fe+3 + 2 H2O = Fe(OH)2+ + 2 H+ - -log_k -5.67 - -delta_h 17.1 kcal - -gamma 5.4 0 -Fe+3 + 3 H2O = Fe(OH)3 + 3 H+ - -log_k -12.56 - -delta_h 24.8 kcal -Fe+3 + 4 H2O = Fe(OH)4- + 4 H+ - -log_k -21.6 - -delta_h 31.9 kcal - -gamma 5.4 0 -Fe+2 + 2H2O = Fe(OH)2 + 2H+ - -log_k -20.57 - -delta_h 28.565 kcal -2 Fe+3 + 2 H2O = Fe2(OH)2+4 + 2 H+ - -log_k -2.95 - -delta_h 13.5 kcal -3 Fe+3 + 4 H2O = Fe3(OH)4+5 + 4 H+ - -log_k -6.3 - -delta_h 14.3 kcal -Fe+3 + Cl- = FeCl+2 - -log_k 1.48 - -delta_h 5.6 kcal - -gamma 5.0 0 -Fe+3 + 2 Cl- = FeCl2+ - -log_k 2.13 - -gamma 5.0 0 -Fe+3 + 3 Cl- = FeCl3 - -log_k 1.13 -Fe+3 + SO4-2 = FeSO4+ - -log_k 4.04 - -delta_h 3.91 kcal - -gamma 5.0 0 -Fe+3 + HSO4- = FeHSO4+2 - -log_k 2.48 -Fe+3 + 2 SO4-2 = Fe(SO4)2- - -log_k 5.38 - -delta_h 4.60 kcal -Fe+3 + HPO4-2 = FeHPO4+ - -log_k 5.43 - -delta_h 5.76 kcal - -gamma 5.0 0 -Fe+3 + H2PO4- = FeH2PO4+2 - -log_k 5.43 - -gamma 5.4 0 -Fe+3 + F- = FeF+2 - -log_k 6.2 - -delta_h 2.7 kcal - -gamma 5.0 0 -Fe+3 + 2 F- = FeF2+ - -log_k 10.8 - -delta_h 4.8 kcal - -gamma 5.0 0 -Fe+3 + 3 F- = FeF3 - -log_k 14.0 - -delta_h 5.4 kcal -Mn+2 + H2O = MnOH+ + H+ - -log_k -10.59 - -delta_h 14.40 kcal - -gamma 5.0 0 -Mn+2 + 3H2O = Mn(OH)3- + 3H+ - -log_k -34.8 - -gamma 5.0 0 -Mn+2 + Cl- = MnCl+ - -log_k 0.61 - -gamma 5.0 0 - -Vm 7.25 -1.08 -25.8 -2.73 3.99 5 0 0 0 1 -Mn+2 + 2 Cl- = MnCl2 - -log_k 0.25 - -Vm 1e-5 0 144 -Mn+2 + 3 Cl- = MnCl3- - -log_k -0.31 - -gamma 5.0 0 - -Vm 11.8 0 0 0 2.4 0 0 0 3.6e-2 1 -Mn+2 + CO3-2 = MnCO3 - -log_k 4.9 -Mn+2 + HCO3- = MnHCO3+ - -log_k 1.95 - -gamma 5.0 0 -Mn+2 + SO4-2 = MnSO4 - -log_k 2.25 - -delta_h 3.370 kcal - -Vm -1.31 -1.83 62.3 -2.7 -Mn+2 + 2 NO3- = Mn(NO3)2 - -log_k 0.6 - -delta_h -0.396 kcal - -Vm 6.16 0 29.4 0 0.9 -Mn+2 + F- = MnF+ - -log_k 0.84 - -gamma 5.0 0 -Mn+2 = Mn+3 + e- - -log_k -25.51 - -delta_h 25.80 kcal - -gamma 9.0 0 -Al+3 + H2O = AlOH+2 + H+ - -log_k -5.0 - -delta_h 11.49 kcal - -analytic -38.253 0.0 -656.27 14.327 - -gamma 5.4 0 - -Vm -1.46 -11.4 10.2 -2.31 1.67 5.4 0 0 0 1 # Barta and Hepler, 1986, Can. J. Chem. 64, 353. -Al+3 + 2 H2O = Al(OH)2+ + 2 H+ - -log_k -10.1 - -delta_h 26.90 kcal - -gamma 5.4 0 - -analytic 88.50 0.0 -9391.6 -27.121 -Al+3 + 3 H2O = Al(OH)3 + 3 H+ - -log_k -16.9 - -delta_h 39.89 kcal - -analytic 226.374 0.0 -18247.8 -73.597 -Al+3 + 4 H2O = Al(OH)4- + 4 H+ - -log_k -22.7 - -delta_h 42.30 kcal - -analytic 51.578 0.0 -11168.9 -14.865 - -gamma 4.5 0 - -dw 1.04e-9 # Mackin & Aller, 1983, GCA 47, 959 -Al+3 + SO4-2 = AlSO4+ - -log_k 3.5 - -delta_h 2.29 kcal - -gamma 4.5 0 -Al+3 + 2SO4-2 = Al(SO4)2- - -log_k 5.0 - -delta_h 3.11 kcal - -gamma 4.5 0 -Al+3 + HSO4- = AlHSO4+2 - -log_k 0.46 -Al+3 + F- = AlF+2 - -log_k 7.0 - -delta_h 1.060 kcal - -gamma 5.4 0 -Al+3 + 2 F- = AlF2+ - -log_k 12.7 - -delta_h 1.980 kcal - -gamma 5.4 0 -Al+3 + 3 F- = AlF3 - -log_k 16.8 - -delta_h 2.160 kcal -Al+3 + 4 F- = AlF4- - -log_k 19.4 - -delta_h 2.20 kcal - -gamma 4.5 0 -# Al+3 + 5 F- = AlF5-2 - # log_k 20.6 - # delta_h 1.840 kcal -# Al+3 + 6 F- = AlF6-3 - # log_k 20.6 - # delta_h -1.670 kcal -H4SiO4 = H3SiO4- + H+ - -log_k -9.83 - -delta_h 6.12 kcal - -analytic -302.3724 -0.050698 15669.69 108.18466 -1119669.0 - -gamma 4 0 - -Vm 7.94 1.0881 5.3224 -2.8240 1.4767 # supcrt + H2O in a1 -H4SiO4 = H2SiO4-2 + 2 H+ - -log_k -23.0 - -delta_h 17.6 kcal - -analytic -294.0184 -0.072650 11204.49 108.18466 -1119669.0 - -gamma 5.4 0 -H4SiO4 + 4 H+ + 6 F- = SiF6-2 + 4 H2O - -log_k 30.18 - -delta_h -16.260 kcal - -gamma 5.0 0 - -Vm 8.5311 13.0492 .6211 -3.3185 2.7716 # supcrt -Ba+2 + H2O = BaOH+ + H+ - -log_k -13.47 - -gamma 5.0 0 -Ba+2 + CO3-2 = BaCO3 - -log_k 2.71 - -delta_h 3.55 kcal - -analytic 0.113 0.008721 - -Vm .2907 -7.0717 8.5295 -2.4867 -.0300 # supcrt -Ba+2 + HCO3- = BaHCO3+ - -log_k 0.982 - -delta_h 5.56 kcal - -analytic -3.0938 0.013669 -Ba+2 + SO4-2 = BaSO4 - -log_k 2.7 -Sr+2 + H2O = SrOH+ + H+ - -log_k -13.29 - -gamma 5.0 0 -Sr+2 + CO3-2 + H+ = SrHCO3+ - -log_k 11.509 - -delta_h 2.489 kcal - -analytic 104.6391 0.04739549 -5151.79 -38.92561 563713.9 - -gamma 5.4 0 -Sr+2 + CO3-2 = SrCO3 - -log_k 2.81 - -delta_h 5.22 kcal - -analytic -1.019 0.012826 - -Vm -.1787 -8.2177 8.9799 -2.4393 -.0300 # supcrt -Sr+2 + SO4-2 = SrSO4 - -log_k 2.29 - -delta_h 2.08 kcal - -Vm 6.7910 -.9666 6.1300 -2.7390 -.0010 # celestite solubility -Li+ + SO4-2 = LiSO4- - -log_k 0.64 - -gamma 5.0 0 -Cu+2 + e- = Cu+ - -log_k 2.72 - -delta_h 1.65 kcal - -gamma 2.5 0 -Cu+ + 2Cl- = CuCl2- - -log_k 5.50 - -delta_h -0.42 kcal - -gamma 4.0 0 -Cu+ + 3Cl- = CuCl3-2 - -log_k 5.70 - -delta_h 0.26 kcal - -gamma 5.0 0.0 -Cu+2 + CO3-2 = CuCO3 - -log_k 6.73 -Cu+2 + 2CO3-2 = Cu(CO3)2-2 - -log_k 9.83 -Cu+2 + HCO3- = CuHCO3+ - -log_k 2.7 -Cu+2 + Cl- = CuCl+ - -log_k 0.43 - -delta_h 8.65 kcal - -gamma 4.0 0 - -Vm -4.19 0 30.4 0 0 4 0 0 1.94e-2 1 -Cu+2 + 2Cl- = CuCl2 - -log_k 0.16 - -delta_h 10.56 kcal - -Vm 26.8 0 -136 -Cu+2 + 3Cl- = CuCl3- - -log_k -2.29 - -delta_h 13.69 kcal - -gamma 4.0 0 -Cu+2 + 4Cl- = CuCl4-2 - -log_k -4.59 - -delta_h 17.78 kcal - -gamma 5.0 0 -Cu+2 + F- = CuF+ - -log_k 1.26 - -delta_h 1.62 kcal -Cu+2 + H2O = CuOH+ + H+ - -log_k -8.0 - -gamma 4.0 0 -Cu+2 + 2 H2O = Cu(OH)2 + 2 H+ - -log_k -13.68 -Cu+2 + 3 H2O = Cu(OH)3- + 3 H+ - -log_k -26.9 -Cu+2 + 4 H2O = Cu(OH)4-2 + 4 H+ - -log_k -39.6 -2Cu+2 + 2H2O = Cu2(OH)2+2 + 2H+ - -log_k -10.359 - -delta_h 17.539 kcal - -analytical 2.497 0.0 -3833.0 -Cu+2 + SO4-2 = CuSO4 - -log_k 2.31 - -delta_h 1.220 kcal - -Vm 5.21 0 -14.6 -Cu+2 + 3HS- = Cu(HS)3- - -log_k 25.9 -Zn+2 + H2O = ZnOH+ + H+ - -log_k -8.96 - -delta_h 13.4 kcal -Zn+2 + 2 H2O = Zn(OH)2 + 2 H+ - -log_k -16.9 -Zn+2 + 3 H2O = Zn(OH)3- + 3 H+ - -log_k -28.4 -Zn+2 + 4 H2O = Zn(OH)4-2 + 4 H+ - -log_k -41.2 -Zn+2 + Cl- = ZnCl+ - -log_k 0.43 - -delta_h 7.79 kcal - -gamma 4.0 0 - -Vm 14.8 -3.91 -105.7 -2.62 0.203 4 0 0 -5.05e-2 1 -Zn+2 + 2 Cl- = ZnCl2 - -log_k 0.45 - -delta_h 8.5 kcal - -Vm -10.1 4.57 241 -2.97 -1e-3 -Zn+2 + 3Cl- = ZnCl3- - -log_k 0.5 - -delta_h 9.56 kcal - -gamma 4.0 0 - -Vm 0.772 15.5 -0.349 -3.42 1.25 0 -7.77 0 0 1 -Zn+2 + 4Cl- = ZnCl4-2 - -log_k 0.2 - -delta_h 10.96 kcal - -gamma 5.0 0 - -Vm 28.42 28 -5.26 -3.94 2.67 0 0 0 4.62e-2 1 -Zn+2 + H2O + Cl- = ZnOHCl + H+ - -log_k -7.48 -Zn+2 + 2HS- = Zn(HS)2 - -log_k 14.94 -Zn+2 + 3HS- = Zn(HS)3- - -log_k 16.1 -Zn+2 + CO3-2 = ZnCO3 - -log_k 5.3 -Zn+2 + 2CO3-2 = Zn(CO3)2-2 - -log_k 9.63 -Zn+2 + HCO3- = ZnHCO3+ - -log_k 2.1 -Zn+2 + SO4-2 = ZnSO4 - -log_k 2.37 - -delta_h 1.36 kcal - -Vm 2.51 0 18.8 -Zn+2 + 2SO4-2 = Zn(SO4)2-2 - -log_k 3.28 - -Vm 10.9 0 -98.7 0 0 0 24 0 -0.236 1 -Zn+2 + Br- = ZnBr+ - -log_k -0.58 -Zn+2 + 2Br- = ZnBr2 - -log_k -0.98 -Zn+2 + F- = ZnF+ - -log_k 1.15 - -delta_h 2.22 kcal -Cd+2 + H2O = CdOH+ + H+ - -log_k -10.08 - -delta_h 13.1 kcal -Cd+2 + 2 H2O = Cd(OH)2 + 2 H+ - -log_k -20.35 -Cd+2 + 3 H2O = Cd(OH)3- + 3 H+ - -log_k -33.3 -Cd+2 + 4 H2O = Cd(OH)4-2 + 4 H+ - -log_k -47.35 -2Cd+2 + H2O = Cd2OH+3 + H+ - -log_k -9.39 - -delta_h 10.9 kcal -Cd+2 + H2O + Cl- = CdOHCl + H+ - -log_k -7.404 - -delta_h 4.355 kcal -Cd+2 + NO3- = CdNO3+ - -log_k 0.4 - -delta_h -5.2 kcal - -Vm 5.95 0 -1.11 0 2.67 7 0 0 1.53e-2 1 -Cd+2 + Cl- = CdCl+ - -log_k 1.98 - -delta_h 0.59 kcal - -Vm 5.69 0 -30.2 0 0 6 0 0 0.112 1 -Cd+2 + 2 Cl- = CdCl2 - -log_k 2.6 - -delta_h 1.24 kcal - -Vm 5.53 -Cd+2 + 3 Cl- = CdCl3- - -log_k 2.4 - -delta_h 3.9 kcal - -Vm 4.6 0 83.9 0 0 0 0 0 0 1 -Cd+2 + CO3-2 = CdCO3 - -log_k 2.9 -Cd+2 + 2CO3-2 = Cd(CO3)2-2 - -log_k 6.4 -Cd+2 + HCO3- = CdHCO3+ - -log_k 1.5 -Cd+2 + SO4-2 = CdSO4 - -log_k 2.46 - -delta_h 1.08 kcal - -Vm 10.4 0 57.9 -Cd+2 + 2SO4-2 = Cd(SO4)2-2 - -log_k 3.5 - -Vm -6.29 0 -93 0 9.5 7 0 0 0 1 -Cd+2 + Br- = CdBr+ - -log_k 2.17 - -delta_h -0.81 kcal -Cd+2 + 2Br- = CdBr2 - -log_k 2.9 -Cd+2 + F- = CdF+ - -log_k 1.1 -Cd+2 + 2F- = CdF2 - -log_k 1.5 -Cd+2 + HS- = CdHS+ - -log_k 10.17 -Cd+2 + 2HS- = Cd(HS)2 - -log_k 16.53 -Cd+2 + 3HS- = Cd(HS)3- - -log_k 18.71 -Cd+2 + 4HS- = Cd(HS)4-2 - -log_k 20.9 -Pb+2 + H2O = PbOH+ + H+ - -log_k -7.71 -Pb+2 + 2 H2O = Pb(OH)2 + 2 H+ - -log_k -17.12 -Pb+2 + 3 H2O = Pb(OH)3- + 3 H+ - -log_k -28.06 -Pb+2 + 4 H2O = Pb(OH)4-2 + 4 H+ - -log_k -39.7 -2 Pb+2 + H2O = Pb2OH+3 + H+ - -log_k -6.36 -Pb+2 + Cl- = PbCl+ - -log_k 1.6 - -delta_h 4.38 kcal - -Vm 2.8934 -.7165 6.0316 -2.7494 .1281 6 # supcrt -Pb+2 + 2 Cl- = PbCl2 - -log_k 1.8 - -delta_h 1.08 kcal - -Vm 6.5402 8.1879 2.5318 -3.1175 -.0300 # supcrt -Pb+2 + 3 Cl- = PbCl3- - -log_k 1.7 - -delta_h 2.17 kcal - -Vm 11.0396 19.1743 -1.7863 -3.5717 .7356 # supcrt -Pb+2 + 4 Cl- = PbCl4-2 - -log_k 1.38 - -delta_h 3.53 kcal - -Vm 16.4150 32.2997 -6.9452 -4.1143 2.3118 # supcrt -Pb+2 + CO3-2 = PbCO3 - -log_k 7.24 -Pb+2 + 2 CO3-2 = Pb(CO3)2-2 - -log_k 10.64 -Pb+2 + HCO3- = PbHCO3+ - -log_k 2.9 -Pb+2 + SO4-2 = PbSO4 - -log_k 2.75 -Pb+2 + 2 SO4-2 = Pb(SO4)2-2 - -log_k 3.47 -Pb+2 + 2HS- = Pb(HS)2 - -log_k 15.27 -Pb+2 + 3HS- = Pb(HS)3- - -log_k 16.57 -3Pb+2 + 4H2O = Pb3(OH)4+2 + 4H+ - -log_k -23.88 - -delta_h 26.5 kcal -Pb+2 + NO3- = PbNO3+ - -log_k 1.17 -Pb+2 + Br- = PbBr+ - -log_k 1.77 - -delta_h 2.88 kcal -Pb+2 + 2Br- = PbBr2 - -log_k 1.44 -Pb+2 + F- = PbF+ - -log_k 1.25 -Pb+2 + 2F- = PbF2 - -log_k 2.56 -Pb+2 + 3F- = PbF3- - -log_k 3.42 -Pb+2 + 4F- = PbF4-2 - -log_k 3.1 - -PHASES -Calcite - CaCO3 = CO3-2 + Ca+2 - -log_k -8.48 - -delta_h -2.297 kcal - -analytic 17.118 -0.046528 -3496 # 0 - 250°C, Ellis, 1959, Plummer and Busenberg, 1982 - -Vm 36.9 cm3/mol # MW (100.09 g/mol) / rho (2.71 g/cm3) -Aragonite - CaCO3 = CO3-2 + Ca+2 - -log_k -8.336 - -delta_h -2.589 kcal - -analytic -171.9773 -0.077993 2903.293 71.595 - -Vm 34.04 -Dolomite - CaMg(CO3)2 = Ca+2 + Mg+2 + 2 CO3-2 - -log_k -17.09 - -delta_h -9.436 kcal - -analytic 31.283 -0.0898 -6438 # 25°C: Hemingway and Robie, 1994; 50–175°C: Bénézeth et al., 2018, GCA 224, 262-275. - -Vm 64.5 -Siderite - FeCO3 = Fe+2 + CO3-2 - -log_k -10.89 - -delta_h -2.480 kcal - -Vm 29.2 -Rhodochrosite - MnCO3 = Mn+2 + CO3-2 - -log_k -11.13 - -delta_h -1.430 kcal - -Vm 31.1 -Strontianite - SrCO3 = Sr+2 + CO3-2 - -log_k -9.271 - -delta_h -0.400 kcal - -analytic 155.0305 0.0 -7239.594 -56.58638 - -Vm 39.69 -Witherite - BaCO3 = Ba+2 + CO3-2 - -log_k -8.562 - -delta_h 0.703 kcal - -analytic 607.642 0.121098 -20011.25 -236.4948 - -Vm 46 -Gypsum - CaSO4:2H2O = Ca+2 + SO4-2 + 2 H2O - -log_k -4.58 - -delta_h -0.109 kcal - -analytic 68.2401 0.0 -3221.51 -25.0627 - -analytical_expression 93.7 5.99E-03 -4e3 -35.019 # better fits the appendix data of Appelo, 2015, AG 55, 62 - -Vm 73.9 # 172.18 / 2.33 (Vm H2O = 13.9 cm3/mol) -Anhydrite - CaSO4 = Ca+2 + SO4-2 - -log_k -4.36 - -delta_h -1.710 kcal - -analytic 84.90 0 -3135.12 -31.79 # 50 - 160oC, 1 - 1e3 atm, anhydrite dissolution, Blount and Dickson, 1973, Am. Mineral. 58, 323. - -Vm 46.1 # 136.14 / 2.95 -Celestite - SrSO4 = Sr+2 + SO4-2 - -log_k -6.63 - -delta_h -4.037 kcal -# -analytic -14805.9622 -2.4660924 756968.533 5436.3588 -40553604.0 - -analytic -7.14 6.11e-3 75 0 0 -1.79e-5 # Howell et al., 1992, JCED 37, 464. - -Vm 46.4 -Barite - BaSO4 = Ba+2 + SO4-2 - -log_k -9.97 - -delta_h 6.35 kcal - -analytical_expression -282.43 -8.972e-2 5822 113.08 # Blount 1977; Templeton, 1960 - -Vm 52.9 -Arcanite - K2SO4 = SO4-2 + 2 K+ - log_k -1.776; -delta_h 5 kcal - -analytical_expression 674.142 0.30423 -18037 -280.236 0 -1.44055e-4 # ref. 3 - # Note, the Linke and Seidell data may give subsaturation in other xpt's, SI = -0.06 - -Vm 65.5 -Mirabilite - Na2SO4:10H2O = SO4-2 + 2 Na+ + 10 H2O - -analytical_expression -301.9326 -0.16232 0 141.078 # ref. 3 - Vm 216 -Thenardite - Na2SO4 = 2 Na+ + SO4-2 - -analytical_expression 57.185 8.6024e-2 0 -30.8341 0 -7.6905e-5 # ref. 3 - -Vm 52.9 -Epsomite - MgSO4:7H2O = Mg+2 + SO4-2 + 7 H2O - log_k -1.74; -delta_h 10.57 kJ - -analytical_expression -3.59 6.21e-3 - Vm 147 -Hexahydrite - MgSO4:6H2O = Mg+2 + SO4-2 + 6 H2O - log_k -1.57; -delta_h 2.35 kJ - -analytical_expression -1.978 1.38e-3 - Vm 132 -Kieserite - MgSO4:H2O = Mg+2 + SO4-2 + H2O - log_k -1.16; -delta_h 9.22 kJ - -analytical_expression 29.485 -5.07e-2 0 -2.662 -7.95e5 - Vm 53.8 -Hydroxyapatite - Ca5(PO4)3OH + 4 H+ = H2O + 3 HPO4-2 + 5 Ca+2 - -log_k -3.421 - -delta_h -36.155 kcal - -Vm 128.9 -Fluorite - CaF2 = Ca+2 + 2 F- - -log_k -10.6 - -delta_h 4.69 kcal - -analytic 66.348 0.0 -4298.2 -25.271 - -Vm 15.7 -SiO2(a) - SiO2 + 2 H2O = H4SiO4 - -log_k -2.71 - -delta_h 3.340 kcal - -analytic -0.26 0.0 -731.0 -Chalcedony - SiO2 + 2 H2O = H4SiO4 - -log_k -3.55 - -delta_h 4.720 kcal - -analytic -0.09 0.0 -1032.0 - -Vm 23.1 -Quartz - SiO2 + 2 H2O = H4SiO4 - -log_k -3.98 - -delta_h 5.990 kcal - -analytic 0.41 0.0 -1309.0 - -Vm 22.67 -Gibbsite - Al(OH)3 + 3 H+ = Al+3 + 3 H2O - -log_k 8.11 - -delta_h -22.800 kcal - -Vm 32.22 -Al(OH)3(a) - Al(OH)3 + 3 H+ = Al+3 + 3 H2O - -log_k 10.8 - -delta_h -26.500 kcal -Kaolinite - Al2Si2O5(OH)4 + 6 H+ = H2O + 2 H4SiO4 + 2 Al+3 - -log_k 7.435 - -delta_h -35.300 kcal - -Vm 99.35 -Albite - NaAlSi3O8 + 8 H2O = Na+ + Al(OH)4- + 3 H4SiO4 - -log_k -18.002 - -delta_h 25.896 kcal - -Vm 101.31 -Anorthite - CaAl2Si2O8 + 8 H2O = Ca+2 + 2 Al(OH)4- + 2 H4SiO4 - -log_k -19.714 - -delta_h 11.580 kcal - -Vm 105.05 -K-feldspar - KAlSi3O8 + 8 H2O = K+ + Al(OH)4- + 3 H4SiO4 - -log_k -20.573 - -delta_h 30.820 kcal - -Vm 108.15 -K-mica - KAl3Si3O10(OH)2 + 10 H+ = K+ + 3 Al+3 + 3 H4SiO4 - -log_k 12.703 - -delta_h -59.376 kcal -Chlorite(14A) - Mg5Al2Si3O10(OH)8 + 16H+ = 5Mg+2 + 2Al+3 + 3H4SiO4 + 6H2O - -log_k 68.38 - -delta_h -151.494 kcal -Ca-Montmorillonite - Ca0.165Al2.33Si3.67O10(OH)2 + 12 H2O = 0.165Ca+2 + 2.33 Al(OH)4- + 3.67 H4SiO4 + 2 H+ - -log_k -45.027 - -delta_h 58.373 kcal - -Vm 156.16 -Talc - Mg3Si4O10(OH)2 + 4 H2O + 6 H+ = 3 Mg+2 + 4 H4SiO4 - -log_k 21.399 - -delta_h -46.352 kcal - -Vm 68.34 -Illite - K0.6Mg0.25Al2.3Si3.5O10(OH)2 + 11.2H2O = 0.6K+ + 0.25Mg+2 + 2.3Al(OH)4- + 3.5H4SiO4 + 1.2H+ - -log_k -40.267 - -delta_h 54.684 kcal - -Vm 141.48 -Chrysotile - Mg3Si2O5(OH)4 + 6 H+ = H2O + 2 H4SiO4 + 3 Mg+2 - -log_k 32.2 - -delta_h -46.800 kcal - -analytic 13.248 0.0 10217.1 -6.1894 - -Vm 106.5808 # 277.11/2.60 -Sepiolite - Mg2Si3O7.5OH:3H2O + 4 H+ + 0.5H2O = 2 Mg+2 + 3 H4SiO4 - -log_k 15.760 - -delta_h -10.700 kcal - -Vm 143.765 -Sepiolite(d) - Mg2Si3O7.5OH:3H2O + 4 H+ + 0.5H2O = 2 Mg+2 + 3 H4SiO4 - -log_k 18.66 -Hematite - Fe2O3 + 6 H+ = 2 Fe+3 + 3 H2O - -log_k -4.008 - -delta_h -30.845 kcal - -Vm 30.39 -Goethite - FeOOH + 3 H+ = Fe+3 + 2 H2O - -log_k -1.0 - -delta_h -14.48 kcal - -Vm 20.84 -Fe(OH)3(a) - Fe(OH)3 + 3 H+ = Fe+3 + 3 H2O - -log_k 4.891 -Pyrite - FeS2 + 2 H+ + 2 e- = Fe+2 + 2 HS- - -log_k -18.479 - -delta_h 11.300 kcal - -Vm 23.48 -FeS(ppt) - FeS + H+ = Fe+2 + HS- - -log_k -3.915 -Mackinawite - FeS + H+ = Fe+2 + HS- - -log_k -4.648 - -Vm 20.45 -Sulfur - S + 2H+ + 2e- = H2S - -log_k 4.882 - -delta_h -9.5 kcal -Vivianite - Fe3(PO4)2:8H2O = 3 Fe+2 + 2 PO4-3 + 8 H2O - -log_k -36.0 -Pyrolusite # H2O added for surface calc's - MnO2:H2O + 4 H+ + 2 e- = Mn+2 + 3 H2O - -log_k 41.38 - -delta_h -65.110 kcal -Hausmannite - Mn3O4 + 8 H+ + 2 e- = 3 Mn+2 + 4 H2O - -log_k 61.03 - -delta_h -100.640 kcal -Manganite - MnOOH + 3 H+ + e- = Mn+2 + 2 H2O - -log_k 25.34 -Pyrochroite - Mn(OH)2 + 2 H+ = Mn+2 + 2 H2O - -log_k 15.2 -Halite - NaCl = Cl- + Na+ - log_k 1.570 - -delta_h 1.37 - #-analytic -713.4616 -.1201241 37302.21 262.4583 -2106915. - -Vm 27.1 -Sylvite - KCl = K+ + Cl- - log_k 0.900 - -delta_h 8.5 - # -analytic 3.984 0.0 -919.55 - Vm 37.5 -# Gases... -CO2(g) - CO2 = CO2 - -log_k -1.468 - -delta_h -4.776 kcal - -analytic 10.5624 -2.3547e-2 -3972.8 0 5.8746e5 1.9194e-5 - -T_c 304.2 # critical T, K - -P_c 72.86 # critical P, atm - -Omega 0.225 # acentric factor -H2O(g) - H2O = H2O - -log_k 1.506; delta_h -44.03 kJ - -T_c 647.3; -P_c 217.60; -Omega 0.344 - -analytic -16.5066 -2.0013E-3 2710.7 3.7646 0 2.24E-6 -O2(g) - O2 = O2 - -log_k -2.8983 - -analytic -7.5001 7.8981e-3 0.0 0.0 2.0027e5 - -T_c 154.6; -P_c 49.80; -Omega 0.021 -H2(g) - H2 = H2 - -log_k -3.1050 - -delta_h -4.184 kJ - -analytic -9.3114 4.6473e-3 -49.335 1.4341 1.2815e5 - -T_c 33.2; -P_c 12.80; -Omega -0.225 -N2(g) - N2 = N2 - -log_k -3.1864 - -analytic -58.453 1.818e-3 3199 17.909 -27460 - -T_c 126.2; -P_c 33.50; -Omega 0.039 -H2S(g) - H2S = H+ + HS- - log_k -7.93 - -delta_h 9.1 - -analytic -45.07 -0.02418 0 17.9205 # H2S solubilities, 0 - 300°C, 1 - 987 atm, Jiang et al., 2020, CG 555, 119816 - -T_c 373.2; -P_c 88.20; -Omega 0.1 -CH4(g) - CH4 = CH4 - -log_k -2.8 - -analytic 10.44 -7.65e-3 -6669 0 1.014e6 # CH4 solubilities 25 - 100°C - -T_c 190.6 ; -P_c 45.40 ; -Omega 0.008 -#Amm(g) -# Amm = Amm -NH3(g) - NH3 = NH3 - -log_k 1.7966 - -analytic -18.758 3.3670e-4 2.5113e3 4.8619 39.192 - -T_c 405.6; -P_c 111.3; -Omega 0.25 -# redox-uncoupled gases -Oxg(g) - Oxg = Oxg - -analytic -7.5001 7.8981e-3 0.0 0.0 2.0027e5 - -T_c 154.6 ; -P_c 49.80 ; -Omega 0.021 -Hdg(g) - Hdg = Hdg - -analytic -9.3114 4.6473e-3 -49.335 1.4341 1.2815e5 - -T_c 33.2 ; -P_c 12.80 ; -Omega -0.225 -Ntg(g) - Ntg = Ntg - -analytic -58.453 1.81800e-3 3199 17.909 -27460 - T_c 126.2 ; -P_c 33.50 ; -Omega 0.039 -Mtg(g) - Mtg = Mtg - -log_k -2.8 - -analytic 10.44 -7.65e-3 -6669 0 1.014e6 # CH4 solubilities 25 - 100°C - -T_c 190.6 ; -P_c 45.40 ; -Omega 0.008 -H2Sg(g) - H2Sg = H+ + HSg- - log_k -7.93 - -delta_h 9.1 - -analytic -45.07 -0.02418 0 17.9205 # H2S solubilities, 0 - 300°C, 1 - 987 atm, Jiang et al., 2020, CG 555, 119816 - -T_c 373.2 ; -P_c 88.20 ; -Omega 0.1 -Melanterite - FeSO4:7H2O = 7 H2O + Fe+2 + SO4-2 - -log_k -2.209 - -delta_h 4.910 kcal - -analytic 1.447 -0.004153 0.0 0.0 -214949.0 -Alunite - KAl3(SO4)2(OH)6 + 6 H+ = K+ + 3 Al+3 + 2 SO4-2 + 6H2O - -log_k -1.4 - -delta_h -50.250 kcal -Jarosite-K - KFe3(SO4)2(OH)6 + 6 H+ = 3 Fe+3 + 6 H2O + K+ + 2 SO4-2 - -log_k -9.21 - -delta_h -31.280 kcal -Zn(OH)2(e) - Zn(OH)2 + 2 H+ = Zn+2 + 2 H2O - -log_k 11.5 -Smithsonite - ZnCO3 = Zn+2 + CO3-2 - -log_k -10.0 - -delta_h -4.36 kcal -Sphalerite - ZnS + H+ = Zn+2 + HS- - -log_k -11.618 - -delta_h 8.250 kcal -Willemite 289 - Zn2SiO4 + 4H+ = 2Zn+2 + H4SiO4 - -log_k 15.33 - -delta_h -33.37 kcal -Cd(OH)2 - Cd(OH)2 + 2 H+ = Cd+2 + 2 H2O - -log_k 13.65 -Otavite 315 - CdCO3 = Cd+2 + CO3-2 - -log_k -12.1 - -delta_h -0.019 kcal -CdSiO3 328 - CdSiO3 + H2O + 2H+ = Cd+2 + H4SiO4 - -log_k 9.06 - -delta_h -16.63 kcal -CdSO4 329 - CdSO4 = Cd+2 + SO4-2 - -log_k -0.1 - -delta_h -14.74 kcal -Cerussite 365 - PbCO3 = Pb+2 + CO3-2 - -log_k -13.13 - -delta_h 4.86 kcal -Anglesite 384 - PbSO4 = Pb+2 + SO4-2 - -log_k -7.79 - -delta_h 2.15 kcal -Pb(OH)2 389 - Pb(OH)2 + 2H+ = Pb+2 + 2H2O - -log_k 8.15 - -delta_h -13.99 kcal - -EXCHANGE_MASTER_SPECIES - X X- -EXCHANGE_SPECIES - X- = X- - -log_k 0.0 - - Na+ + X- = NaX - -log_k 0.0 - -gamma 4.08 0.082 - - K+ + X- = KX - -log_k 0.7 - -gamma 3.5 0.015 - -delta_h -4.3 # Jardine & Sparks, 1984 - - Li+ + X- = LiX - -log_k -0.08 - -gamma 6.0 0 - -delta_h 1.4 # Merriam & Thomas, 1956 - -# !!!!! -# H+ + X- = HX -# -log_k 1.0 -# -gamma 9.0 0 - -# AmmH+ + X- = AmmHX - NH4+ + X- = NH4X - -log_k 0.6 - -gamma 2.5 0 - -delta_h -2.4 # Laudelout et al., 1968 - - Ca+2 + 2X- = CaX2 - -log_k 0.8 - -gamma 5.0 0.165 - -delta_h 7.2 # Van Bladel & Gheyl, 1980 - - Mg+2 + 2X- = MgX2 - -log_k 0.6 - -gamma 5.5 0.2 - -delta_h 7.4 # Laudelout et al., 1968 - - Sr+2 + 2X- = SrX2 - -log_k 0.91 - -gamma 5.26 0.121 - -delta_h 5.5 # Laudelout et al., 1968 - - Ba+2 + 2X- = BaX2 - -log_k 0.91 - -gamma 4.0 0.153 - -delta_h 4.5 # Laudelout et al., 1968 - - Mn+2 + 2X- = MnX2 - -log_k 0.52 - -gamma 6.0 0 - - Fe+2 + 2X- = FeX2 - -log_k 0.44 - -gamma 6.0 0 - - Cu+2 + 2X- = CuX2 - -log_k 0.6 - -gamma 6.0 0 - - Zn+2 + 2X- = ZnX2 - -log_k 0.8 - -gamma 5.0 0 - - Cd+2 + 2X- = CdX2 - -log_k 0.8 - -gamma 0.0 0 - - Pb+2 + 2X- = PbX2 - -log_k 1.05 - -gamma 0.0 0 - - Al+3 + 3X- = AlX3 - -log_k 0.41 - -gamma 9.0 0 - - AlOH+2 + 2X- = AlOHX2 - -log_k 0.89 - -gamma 0.0 0 - -SURFACE_MASTER_SPECIES - Hfo_s Hfo_sOH - Hfo_w Hfo_wOH -SURFACE_SPECIES -# All surface data from -# Dzombak and Morel, 1990 -# -# -# Acid-base data from table 5.7 -# -# strong binding site--Hfo_s, - - Hfo_sOH = Hfo_sOH - -log_k 0 - - Hfo_sOH + H+ = Hfo_sOH2+ - -log_k 7.29 # = pKa1,int - - Hfo_sOH = Hfo_sO- + H+ - -log_k -8.93 # = -pKa2,int - -# weak binding site--Hfo_w - - Hfo_wOH = Hfo_wOH - -log_k 0 - - Hfo_wOH + H+ = Hfo_wOH2+ - -log_k 7.29 # = pKa1,int - - Hfo_wOH = Hfo_wO- + H+ - -log_k -8.93 # = -pKa2,int -############################################### -# CATIONS # -############################################### -# -# Cations from table 10.1 or 10.5 -# -# Calcium - Hfo_sOH + Ca+2 = Hfo_sOHCa+2 - -log_k 4.97 - - Hfo_wOH + Ca+2 = Hfo_wOCa+ + H+ - -log_k -5.85 -# Strontium - Hfo_sOH + Sr+2 = Hfo_sOHSr+2 - -log_k 5.01 - - Hfo_wOH + Sr+2 = Hfo_wOSr+ + H+ - -log_k -6.58 - - Hfo_wOH + Sr+2 + H2O = Hfo_wOSrOH + 2H+ - -log_k -17.6 -# Barium - Hfo_sOH + Ba+2 = Hfo_sOHBa+2 - -log_k 5.46 - - Hfo_wOH + Ba+2 = Hfo_wOBa+ + H+ - -log_k -7.2 # table 10.5 -# -# Cations from table 10.2 -# -# Cadmium - Hfo_sOH + Cd+2 = Hfo_sOCd+ + H+ - -log_k 0.47 - - Hfo_wOH + Cd+2 = Hfo_wOCd+ + H+ - -log_k -2.91 -# Zinc - Hfo_sOH + Zn+2 = Hfo_sOZn+ + H+ - -log_k 0.99 - - Hfo_wOH + Zn+2 = Hfo_wOZn+ + H+ - -log_k -1.99 -# Copper - Hfo_sOH + Cu+2 = Hfo_sOCu+ + H+ - -log_k 2.89 - - Hfo_wOH + Cu+2 = Hfo_wOCu+ + H+ - -log_k 0.6 # table 10.5 -# Lead - Hfo_sOH + Pb+2 = Hfo_sOPb+ + H+ - -log_k 4.65 - - Hfo_wOH + Pb+2 = Hfo_wOPb+ + H+ - -log_k 0.3 # table 10.5 -# -# Derived constants table 10.5 -# -# Magnesium - Hfo_wOH + Mg+2 = Hfo_wOMg+ + H+ - -log_k -4.6 -# Manganese - Hfo_sOH + Mn+2 = Hfo_sOMn+ + H+ - -log_k -0.4 # table 10.5 - - Hfo_wOH + Mn+2 = Hfo_wOMn+ + H+ - -log_k -3.5 # table 10.5 -# Iron, strong site: Appelo, Van der Weiden, Tournassat & Charlet, EST 36, 3096 - Hfo_sOH + Fe+2 = Hfo_sOFe+ + H+ - -log_k -0.95 -# Iron, weak site: Liger et al., GCA 63, 2939, re-optimized for D&M - Hfo_wOH + Fe+2 = Hfo_wOFe+ + H+ - -log_k -2.98 - - Hfo_wOH + Fe+2 + H2O = Hfo_wOFeOH + 2H+ - -log_k -11.55 -############################################### -# ANIONS # -############################################### -# -# Anions from table 10.6 -# -# Phosphate - Hfo_wOH + PO4-3 + 3H+ = Hfo_wH2PO4 + H2O - -log_k 31.29 - - Hfo_wOH + PO4-3 + 2H+ = Hfo_wHPO4- + H2O - -log_k 25.39 - - Hfo_wOH + PO4-3 + H+ = Hfo_wPO4-2 + H2O - -log_k 17.72 -# -# Anions from table 10.7 -# -# Borate - Hfo_wOH + H3BO3 = Hfo_wH2BO3 + H2O - -log_k 0.62 -# -# Anions from table 10.8 -# -# Sulfate - Hfo_wOH + SO4-2 + H+ = Hfo_wSO4- + H2O - -log_k 7.78 - - Hfo_wOH + SO4-2 = Hfo_wOHSO4-2 - -log_k 0.79 -# -# Derived constants table 10.10 -# - Hfo_wOH + F- + H+ = Hfo_wF + H2O - -log_k 8.7 - - Hfo_wOH + F- = Hfo_wOHF- - -log_k 1.6 -# -# Carbonate: Van Geen et al., 1994 reoptimized for D&M model -# - Hfo_wOH + CO3-2 + H+ = Hfo_wCO3- + H2O - -log_k 12.56 - - Hfo_wOH + CO3-2 + 2H+= Hfo_wHCO3 + H2O - -log_k 20.62 -# -# Silicate: Swedlund, P.J. and Webster, J.G., 1999. Water Research 33, 3413-3422. -# - Hfo_wOH + H4SiO4 = Hfo_wH3SiO4 + H2O ; log_K 4.28 - Hfo_wOH + H4SiO4 = Hfo_wH2SiO4- + H+ + H2O ; log_K -3.22 - Hfo_wOH + H4SiO4 = Hfo_wHSiO4-2 + 2H+ + H2O ; log_K -11.69 - -RATES - -########### -#Quartz -########### -# -####### -# Example of quartz kinetic rates block: -# KINETICS -# Quartz -# -m0 158.8 # 90 % Qu -# -parms 0.146 1.5 -# -step 3.1536e8 in 10 -# -tol 1e-12 - -Quartz - -start -1 REM Specific rate k from Rimstidt and Barnes, 1980, GCA 44,1683 -2 REM k = 10^-13.7 mol/m2/s (25 C), Ea = 90 kJ/mol -3 REM sp. rate * parm(2) due to salts (Dove and Rimstidt, MSA Rev. 29, 259) -4 REM PARM(1) = Specific area of Quartz, m^2/mol Quartz -5 REM PARM(2) = salt correction: (1 + 1.5 * c_Na (mM)), < 35 - -10 dif_temp = 1/TK - 1/298 -20 pk_w = 13.7 + 4700.4 * dif_temp -40 moles = PARM(1) * M0 * PARM(2) * (M/M0)^0.67 * 10^-pk_w * (1 - SR("Quartz")) -# Integrate... -50 SAVE moles * TIME - -end - -########### -#K-feldspar -########### -# -# Sverdrup and Warfvinge, 1995, Estimating field weathering rates -# using laboratory kinetics: Reviews in mineralogy and geochemistry, -# vol. 31, p. 485-541. -# -# As described in: -# Appelo and Postma, 2005, Geochemistry, groundwater -# and pollution, 2nd Edition: A.A. Balkema Publishers, -# p. 162-163 and 395-399. -# -# Assume soil is 10% K-feldspar by mass in 1 mm spheres (radius 0.05 mm) -# Assume density of rock and Kspar is 2600 kg/m^3 = 2.6 kg/L -# GFW Kspar 0.278 kg/mol -# -# Moles of Kspar per liter pore space calculation: -# Mass of rock per liter pore space = 0.7*2.6/0.3 = 6.07 kg rock/L pore space -# Mass of Kspar per liter pore space 6.07x0.1 = 0.607 kg Kspar/L pore space -# Moles of Kspar per liter pore space 0.607/0.278 = 2.18 mol Kspar/L pore space -# -# Specific area calculation: -# Volume of sphere 4/3 x pi x r^3 = 5.24e-13 m^3 Kspar/sphere -# Mass of sphere 2600 x 5.24e-13 = 1.36e-9 kg Kspar/sphere -# Moles of Kspar in sphere 1.36e-9/0.278 = 4.90e-9 mol Kspar/sphere -# Surface area of one sphere 4 x pi x r^2 = 3.14e-8 m^2/sphere -# Specific area of K-feldspar in sphere 3.14e-8/4.90e-9 = 6.41 m^2/mol Kspar -# -# -# Example of KINETICS data block for K-feldspar rate: -# KINETICS 1 -# K-feldspar -# -m0 2.18 # 10% Kspar, 0.1 mm cubes -# -m 2.18 # Moles per L pore space -# -parms 6.41 0.1 # m^2/mol Kspar, fraction adjusts lab rate to field rate -# -time 1.5 year in 40 - -K-feldspar - -start -1 REM Sverdrup and Warfvinge, 1995, mol m^-2 s^-1 -2 REM PARM(1) = Specific area of Kspar m^2/mol Kspar -3 REM PARM(2) = Adjusts lab rate to field rate -4 REM temp corr: from A&P, p. 162. E (kJ/mol) / R / 2.303 = H in H*(1/T-1/281) -5 REM K-Feldspar parameters -10 DATA 11.7, 0.5, 4e-6, 0.4, 500e-6, 0.15, 14.5, 0.14, 0.15, 13.1, 0.3 -20 RESTORE 10 -30 READ pK_H, n_H, lim_Al, x_Al, lim_BC, x_BC, pK_H2O, z_Al, z_BC, pK_OH, o_OH -40 DATA 3500, 2000, 2500, 2000 -50 RESTORE 40 -60 READ e_H, e_H2O, e_OH, e_CO2 -70 pk_CO2 = 13 -80 n_CO2 = 0.6 -100 REM Generic rate follows -110 dif_temp = 1/TK - 1/281 -120 BC = ACT("Na+") + ACT("K+") + ACT("Mg+2") + ACT("Ca+2") -130 REM rate by H+ -140 pk_H = pk_H + e_H * dif_temp -150 rate_H = 10^-pk_H * ACT("H+")^n_H / ((1 + ACT("Al+3") / lim_Al)^x_Al * (1 + BC / lim_BC)^x_BC) -160 REM rate by hydrolysis -170 pk_H2O = pk_H2O + e_H2O * dif_temp -180 rate_H2O = 10^-pk_H2O / ((1 + ACT("Al+3") / lim_Al)^z_Al * (1 + BC / lim_BC)^z_BC) -190 REM rate by OH- -200 pk_OH = pk_OH + e_OH * dif_temp -210 rate_OH = 10^-pk_OH * ACT("OH-")^o_OH -220 REM rate by CO2 -230 pk_CO2 = pk_CO2 + e_CO2 * dif_temp -240 rate_CO2 = 10^-pk_CO2 * (SR("CO2(g)"))^n_CO2 -250 rate = rate_H + rate_H2O + rate_OH + rate_CO2 -260 area = PARM(1) * M0 *(M/M0)^0.67 -270 rate = PARM(2) * area * rate * (1-SR("K-feldspar")) -280 moles = rate * TIME -290 SAVE moles - -end - - -########### -#Albite -########### -# -# Sverdrup and Warfvinge, 1995, Estimating field weathering rates -# using laboratory kinetics: Reviews in mineralogy and geochemistry, -# vol. 31, p. 485-541. -# -# As described in: -# Appelo and Postma, 2005, Geochemistry, groundwater -# and pollution, 2nd Edition: A.A. Balkema Publishers, -# p. 162-163 and 395-399. -# -# Example of KINETICS data block for Albite rate: -# KINETICS 1 -# Albite -# -m0 0.46 # 2% Albite, 0.1 mm cubes -# -m 0.46 # Moles per L pore space -# -parms 6.04 0.1 # m^2/mol Albite, fraction adjusts lab rate to field rate -# -time 1.5 year in 40 -# -# Assume soil is 2% Albite by mass in 1 mm spheres (radius 0.05 mm) -# Assume density of rock and Albite is 2600 kg/m^3 = 2.6 kg/L -# GFW Albite 0.262 kg/mol -# -# Moles of Albite per liter pore space calculation: -# Mass of rock per liter pore space = 0.7*2.6/0.3 = 6.07 kg rock/L pore space -# Mass of Albite per liter pore space 6.07x0.02 = 0.121 kg Albite/L pore space -# Moles of Albite per liter pore space 0.607/0.262 = 0.46 mol Albite/L pore space -# -# Specific area calculation: -# Volume of sphere 4/3 x pi x r^3 = 5.24e-13 m^3 Albite/sphere -# Mass of sphere 2600 x 5.24e-13 = 1.36e-9 kg Albite/sphere -# Moles of Albite in sphere 1.36e-9/0.262 = 5.20e-9 mol Albite/sphere -# Surface area of one sphere 4 x pi x r^2 = 3.14e-8 m^2/sphere -# Specific area of Albite in sphere 3.14e-8/5.20e-9 = 6.04 m^2/mol Albite - -Albite - -start -1 REM Sverdrup and Warfvinge, 1995, mol m^-2 s^-1 -2 REM PARM(1) = Specific area of Albite m^2/mol Albite -3 REM PARM(2) = Adjusts lab rate to field rate -4 REM temp corr: from A&P, p. 162. E (kJ/mol) / R / 2.303 = H in H*(1/T-1/281) -5 REM Albite parameters -10 DATA 11.5, 0.5, 4e-6, 0.4, 500e-6, 0.2, 13.7, 0.14, 0.15, 11.8, 0.3 -20 RESTORE 10 -30 READ pK_H, n_H, lim_Al, x_Al, lim_BC, x_BC, pK_H2O, z_Al, z_BC, pK_OH, o_OH -40 DATA 3500, 2000, 2500, 2000 -50 RESTORE 40 -60 READ e_H, e_H2O, e_OH, e_CO2 -70 pk_CO2 = 13 -80 n_CO2 = 0.6 -100 REM Generic rate follows -110 dif_temp = 1/TK - 1/281 -120 BC = ACT("Na+") + ACT("K+") + ACT("Mg+2") + ACT("Ca+2") -130 REM rate by H+ -140 pk_H = pk_H + e_H * dif_temp -150 rate_H = 10^-pk_H * ACT("H+")^n_H / ((1 + ACT("Al+3") / lim_Al)^x_Al * (1 + BC / lim_BC)^x_BC) -160 REM rate by hydrolysis -170 pk_H2O = pk_H2O + e_H2O * dif_temp -180 rate_H2O = 10^-pk_H2O / ((1 + ACT("Al+3") / lim_Al)^z_Al * (1 + BC / lim_BC)^z_BC) -190 REM rate by OH- -200 pk_OH = pk_OH + e_OH * dif_temp -210 rate_OH = 10^-pk_OH * ACT("OH-")^o_OH -220 REM rate by CO2 -230 pk_CO2 = pk_CO2 + e_CO2 * dif_temp -240 rate_CO2 = 10^-pk_CO2 * (SR("CO2(g)"))^n_CO2 -250 rate = rate_H + rate_H2O + rate_OH + rate_CO2 -260 area = PARM(1) * M0 *(M/M0)^0.67 -270 rate = PARM(2) * area * rate * (1-SR("Albite")) -280 moles = rate * TIME -290 SAVE moles - -end - -######## -#Calcite -######## -# Example of KINETICS data block for calcite rate, -# in mmol/cm2/s, Plummer et al., 1978, AJS 278, 179; Appelo et al., AG 13, 257. -# KINETICS 1 -# Calcite -# -tol 1e-8 -# -m0 3.e-3 -# -m 3.e-3 -# -parms 1.67e5 0.6 # cm^2/mol calcite, exp factor -# -time 1 day - -Calcite - -start -1 REM PARM(1) = specific surface area of calcite, cm^2/mol calcite -2 REM PARM(2) = exponent for M/M0 - -10 si_cc = SI("Calcite") -20 IF (M <= 0 and si_cc < 0) THEN GOTO 200 -30 k1 = 10^(0.198 - 444.0 / TK ) -40 k2 = 10^(2.84 - 2177.0 /TK ) -50 IF TC <= 25 THEN k3 = 10^(-5.86 - 317.0 / TK) -60 IF TC > 25 THEN k3 = 10^(-1.1 - 1737.0 / TK ) -80 IF M0 > 0 THEN area = PARM(1)*M0*(M/M0)^PARM(2) ELSE area = PARM(1)*M -110 rate = area * (k1 * ACT("H+") + k2 * ACT("CO2") + k3 * ACT("H2O")) -120 rate = rate * (1 - 10^(2/3*si_cc)) -130 moles = rate * 0.001 * TIME # convert from mmol to mol -200 SAVE moles - -end - -####### -#Pyrite -####### -# -# Williamson, M.A. and Rimstidt, J.D., 1994, -# Geochimica et Cosmochimica Acta, v. 58, p. 5443-5454, -# rate equation is mol m^-2 s^-1. -# -# Example of KINETICS data block for pyrite rate: -# KINETICS 1 -# Pyrite -# -tol 1e-8 -# -m0 5.e-4 -# -m 5.e-4 -# -parms 0.3 0.67 .5 -0.11 -# -time 1 day in 10 -Pyrite - -start -1 REM Williamson and Rimstidt, 1994 -2 REM PARM(1) = log10(specific area), log10(m^2 per mole pyrite) -3 REM PARM(2) = exp for (M/M0) -4 REM PARM(3) = exp for O2 -5 REM PARM(4) = exp for H+ - -10 REM Dissolution in presence of DO -20 if (M <= 0) THEN GOTO 200 -30 if (SI("Pyrite") >= 0) THEN GOTO 200 -40 log_rate = -8.19 + PARM(3)*LM("O2") + PARM(4)*LM("H+") -50 log_area = PARM(1) + LOG10(M0) + PARM(2)*LOG10(M/M0) -60 moles = 10^(log_area + log_rate) * TIME -200 SAVE moles - -end - -########## -#Organic_C -########## -# -# Example of KINETICS data block for SOC (sediment organic carbon): -# KINETICS 1 -# Organic_C -# -formula C -# -tol 1e-8 -# -m 5e-3 # SOC in mol -# -time 30 year in 15 -Organic_C - -start -1 REM Additive Monod kinetics for SOC (sediment organic carbon) -2 REM Electron acceptors: O2, NO3, and SO4 - -10 if (M <= 0) THEN GOTO 200 -20 mO2 = MOL("O2") -30 mNO3 = TOT("N(5)") -40 mSO4 = TOT("S(6)") -50 k_O2 = 1.57e-9 # 1/sec -60 k_NO3 = 1.67e-11 # 1/sec -70 k_SO4 = 1.e-13 # 1/sec -80 rate = k_O2 * mO2/(2.94e-4 + mO2) -90 rate = rate + k_NO3 * mNO3/(1.55e-4 + mNO3) -100 rate = rate + k_SO4 * mSO4/(1.e-4 + mSO4) -110 moles = rate * M * (M/M0) * TIME -200 SAVE moles - -end - -########### -#Pyrolusite -########### -# -# Postma, D. and Appelo, C.A.J., 2000, GCA, vol. 64, pp. 1237-1247. -# Rate equation given as mol L^-1 s^-1 -# -# Example of KINETICS data block for Pyrolusite -# KINETICS 1-12 -# Pyrolusite -# -tol 1.e-7 -# -m0 0.1 -# -m 0.1 -# -time 0.5 day in 10 -Pyrolusite - -start -10 if (M <= 0) THEN GOTO 200 -20 sr_pl = SR("Pyrolusite") -30 if (sr_pl > 1) THEN GOTO 100 -40 REM sr_pl <= 1, undersaturated -50 Fe_t = TOT("Fe(2)") -60 if Fe_t < 1e-8 then goto 200 -70 moles = 6.98e-5 * Fe_t * (M/M0)^0.67 * TIME * (1 - sr_pl) -80 GOTO 200 -100 REM sr_pl > 1, supersaturated -110 moles = 2e-3 * 6.98e-5 * (1 - sr_pl) * TIME -200 SAVE moles * SOLN_VOL - -end -END -# ============================================================================================= -#(a) means amorphous. (d) means disordered, or less crystalline. -#(14A) refers to 14 angstrom spacing of clay planes. FeS(ppt), -#precipitated, indicates an initial precipitate that is less crystalline. -#Zn(OH)2(e) indicates a specific crystal form, epsilon. -# ============================================================================================= -# For the reaction aA + bB = cC + dD, -# with delta_v = c*Vm(C) + d*Vm(D) - a*Vm(A) - b*Vm(B), -# PHREEQC adds the pressure term to log_k: -= delta_v * (P - 1) / (2.3RT). -# Vm(A) is volume of A, cm3/mol, P is pressure, atm, R is the gas constant, T is Kelvin. -# Gas-pressures and fugacity coefficients are calculated with Peng-Robinson's EOS. -# Binary interaction coefficients from Soreide and Whitson, 1992, FPE 77, 217 are -# hard-coded in calc_PR(): -# kij CH4 CO2 H2S N2 -# H2O 0.49 0.19 0.19 0.49 -# ============================================================================================= -# The molar volumes of solids are entered with -# -Vm vm cm3/mol -# vm is the molar volume, cm3/mol (default), but dm3/mol and m3/mol are permitted. -# Data for minerals' vm (= MW (g/mol) / rho (g/cm3)) are defined using rho from -# Deer, Howie and Zussman, The rock-forming minerals, Longman. -# -------------------- -# Temperature- and pressure-dependent volumina of aqueous species are calculated with a Redlich- -# type equation (cf. Redlich and Meyer, Chem. Rev. 64, 221), from parameters entered with -# -Vm a1 a2 a3 a4 W a0 i1 i2 i3 i4 -# The volume (cm3/mol) is -# Vm(T, pb, I) = 41.84 * (a1 * 0.1 + a2 * 100 / (2600 + pb) + a3 / (T - 228) + -# a4 * 1e4 / (2600 + pb) / (T - 228) - W * QBrn) -# + z^2 / 2 * Av * f(I^0.5) -# + (i1 + i2 / (T - 228) + i3 * (T - 228)) * I^i4 -# Volumina at I = 0 are obtained using supcrt92 formulas (Johnson et al., 1992, CG 18, 899). -# 41.84 transforms cal/bar/mol into cm3/mol. -# pb is pressure in bar. -# W * QBrn is the energy of solvation, calculated from W and the pressure dependence of the Born equation, -# W is fitted on measured solution densities. -# z is charge of the solute species. -# Av is the Debye-Hückel limiting slope (DH_AV in PHREEQC basic). -# a0 is the ion-size parameter in the extended Debye-Hückel equation: -# f(I^0.5) = I^0.5 / (1 + a0 * DH_B * I^0.5), -# a0 = -gamma x for cations, = 0 for anions. -# For details, consult ref. 1. -# ============================================================================================= -# The viscosity is calculated with a (modified) Jones-Dole equation: -# viscos / viscos_0 = 1 + A Sum(0.5 z_i m_i) + fan (B_i m_i + D_i m_i n_i) -# Parameters are for calculating the B and D terms: -# -viscosity 9.35e-2 -8.31e-2 2.487e-2 4.49e-4 2.01e-2 1.570 0 -# # b0 b1 b2 d1 d2 d3 tan -# z_i is absolute charge number, m_i is molality of i -# B_i = b0 + b1 exp(-b2 * tc) -# fan = (2 - tan V_i / V_Cl-), corrects for the volume of anions -# D_i = d1 + exp(-d2 tc) -# n_i = ((1 + fI)^d3 + ((z_i^2 + z_i) / 2 · m_i)d^3 / (2 + fI), fI is an ionic strength term. -# For details, consult ref. 4. -# -# ref. 1: Appelo, Parkhurst and Post, 2014. Geochim. Cosmochim. Acta 125, 49–67. -# ref. 2: Procedures from ref. 1 using data compiled by Laliberté, 2009, J. Chem. Eng. Data 54, 1725. -# ref. 3: Appelo, 2017, Cem. Concr. Res. 101, 102-113. -# ref. 4: Appelo and Parkhurst in prep., for details see subroutine viscosity in transport.cpp -# -# ============================================================================================= -# It remains the responsibility of the user to check the calculated results, for example with -# measured solubilities as a function of (P, T). diff --git a/pybind/tests/test_basic.py b/pybind/tests/test_basic.py deleted file mode 100644 index 33990eac1..000000000 --- a/pybind/tests/test_basic.py +++ /dev/null @@ -1,396 +0,0 @@ -import os -##import phreeqcrm as rm -import numpy as np -from numpy.testing import assert_array_almost_equal, assert_array_less -import pytest - -from phreeqcrm import bmi_phreeqcrm, IRM_RESULT - -from constants import FilePaths - -def test_main(): - # for debugging - print(f"PYTHONPATH={os.getenv('PYTHONPATH')}") - -def test_dtor(): - model = bmi_phreeqcrm() - del model - -def test_initialized(): - model = bmi_phreeqcrm() - assert(not model._initialized) - -def test_finalize_not_initialized(): - model = bmi_phreeqcrm() - - model.finalize() - assert(not model._initialized) - -def test_get_components_is_tuple(): - model = bmi_phreeqcrm() - - # may want to make this throw to be consistent - components = model.get_components() - assert(isinstance(components, tuple)) - assert(len(components) == 0) - -def test_get_grid_size_raises_uninitialized(): - model = bmi_phreeqcrm() - with pytest.raises(RuntimeError, match="must call initialize first"): - model.get_grid_size(0) - -def test_get_input_var_names_raises_uninitialized(): - model = bmi_phreeqcrm() - with pytest.raises(RuntimeError, match="must call initialize first"): - model.get_input_var_names() - -def test_get_output_var_names_raises_uninitialized(): - model = bmi_phreeqcrm() - with pytest.raises(RuntimeError, match="must call initialize first"): - model.get_output_var_names() - -def test_get_value_ptr_raises_uninitialized(): - model = bmi_phreeqcrm() - with pytest.raises(RuntimeError, match="must call initialize first"): - model.get_value_ptr("Temperature") - -def test_update_raises_uninitialized(): - model = bmi_phreeqcrm() - with pytest.raises(RuntimeError, match="must call initialize first"): - model.update() - -def test_update_until_raises_uninitialized(): - model = bmi_phreeqcrm() - with pytest.raises(RuntimeError, match="must call initialize first"): - model.update_until(1.0) - -def test_initialize_AdvectBMI(): - model = bmi_phreeqcrm() - - assert(not model._initialized) - model.initialize(FilePaths.YAML) - assert(model._initialized) - -def test_get_grid_size_AdvectBMI(): - model = bmi_phreeqcrm() - model.initialize(FilePaths.YAML) - - nxyz = model.get_grid_size(0) - assert(nxyz == 40) - -def test_get_value_ptr_is_ndarray(): - model = bmi_phreeqcrm() - model.initialize(FilePaths.YAML) - nxyz = model.get_grid_size(0) - - temperature = model.get_value_ptr("Temperature") - assert(isinstance(temperature, np.ndarray)) - assert(len(temperature) == nxyz) - assert(temperature[0] == pytest.approx(25.)) - assert(temperature[1] == pytest.approx(25.)) - assert(temperature[nxyz - 2] == pytest.approx(25.)) - assert(temperature[nxyz - 2] == pytest.approx(25.)) - -def test_get_temperature_ptr_is_writeable(): - model = bmi_phreeqcrm() - model.initialize(FilePaths.YAML) - nxyz = model.get_grid_size(0) - - # test WRITEABLE - temperature = model.get_value_ptr("Temperature") - assert(temperature.flags.writeable) - - # write - temperature[0] = 0.0 - - temperature = model.get_value_ptr("Temperature") - assert(isinstance(temperature, np.ndarray)) - assert(len(temperature) == nxyz) - assert(temperature[0] == pytest.approx(0.0)) - assert(temperature[1] == pytest.approx(25.)) - assert(temperature[nxyz - 2] == pytest.approx(25.)) - assert(temperature[nxyz - 1] == pytest.approx(25.)) - -def test_get_ComponentCount_ptr_is_readonly(): - model = bmi_phreeqcrm() - model.initialize(FilePaths.YAML) - - # make sure get_value_ptr can be called before get_input_var_names - component_count = model.get_value_ptr("ComponentCount") - assert(isinstance(component_count, np.ndarray)) - assert(len(component_count) == 1) - assert(component_count.size == 1) - assert(component_count[0] == 8) - assert(component_count.dtype == 'int32') - - # make sure ComponentCount is readonly - with pytest.raises(ValueError, match="assignment destination is read-only"): - component_count[0] = 25 - assert(component_count[0] == 8) - -def test_get_input_var_names_is_tuple(): - model = bmi_phreeqcrm() - model.initialize(FilePaths.YAML) - - input_vars = model.get_input_var_names() - assert(isinstance(input_vars, tuple)) - assert(isinstance(input_vars[0], str)) - -def test_get_output_var_names_is_tuple(): - model = bmi_phreeqcrm() - model.initialize(FilePaths.YAML) - - output_vars = model.get_output_var_names() - assert(isinstance(output_vars, tuple)) - assert(isinstance(output_vars[0], str)) - -def test_get_components_is_tuple(): - model = bmi_phreeqcrm() - model.initialize(FilePaths.YAML) - - components = model.get_components() - assert(isinstance(components, tuple)) - assert(len(components) == 8) - -def test_get_Components_ptr_is_ndarray(): - model = bmi_phreeqcrm() - model.initialize(FilePaths.YAML) - - components = model.get_value_ptr("Components") - assert(isinstance(components, np.ndarray)) - assert(components.dtype == '