Merge pull request #2 from avdudchenko/main

Recommiting first API

.github/workflows/checks.yml

@@ -9,7 +9,7 @@ on:
 env:
   PYTEST_ADDOPTS: --color=yes
   PIP_PROGRESS_BAR: "off"
-
+  KERNEL_NAME: reaktoro-pse-dev
 defaults:
   run:
     # -l: login shell, needed when using Conda:
@@ -83,17 +83,29 @@ jobs:
         name: Install (dev)
         run: |
           pip install -r requirements-dev.txt
+          pip list
       - if: matrix.install-mode == 'standard'
         name: Install (standard)
         run: |
           pip install "git+${{ format('{0}/{1}@{2}', github.server_url, github.repository, github.ref) }}"
+      - name: Install (idaes-solvers)
+        run: |
+          idaes get-extensions
       - if: matrix.coverage
         name: Enable coverage for pytest
         run: echo PYTEST_ADDOPTS="$PYTEST_ADDOPTS --cov --cov-report=xml" >> $GITHUB_ENV
       - name: Run pytest
         run: |
           pip install pytest  # ensure pytest is installed (should do nothing if already present from requirements-dev.txt)
           pytest --pyargs reaktoro_pse --verbose
+      - name: Install Jupyter kernel
+        run: |
+          jupyter kernelspec list
+          python -m ipykernel install --user --name "${{ env.KERNEL_NAME }}"
+          jupyter kernelspec list
+      - name: Run pytest with nbmake
+        run: |
+          pytest --nbmake --nbmake-kernel="${{ env.KERNEL_NAME }}" src/reaktoro_pse/tutorials/*.ipynb
       - name: Upload coverage report as job artifact
         if: matrix.coverage
         uses: actions/upload-artifact@v4

.gitignore

@@ -160,3 +160,6 @@ cython_debug/
 #  and can be added to the global gitignore or merged into this file.  For a more nuclear
 #  option (not recommended) you can uncomment the following to ignore the entire idea folder.
 #.idea/
+*.jpg
+/.vscode
+src/reaktoro_pse/copyright.txt

README.md

@@ -1,15 +1,117 @@
-# reaktoro-pse
+# reaktoro-pse introduction
+## 1. Overview
+This is a package for configuring [Reaktoro](https://reaktoro.org/index.html) as a gray box model in [Pyomo](https://pyomo.readthedocs.io/en/stable/), [IDAES-PSE](https://idaes-pse.readthedocs.io/en/stable/), and [WaterTAP](https://watertap.readthedocs.io/en/stable/) modeling libraries. This package is not meant to replace or act as a higher level API for Reaktoro - it is only meant to enable setting up Reaktoro equilibrium problems as blocks on Pyomo models and automate transferring Reaktoro data into Pyomo variables. 
 
-## Compatibility
+## 2. Prerequisites
+**The user must familiarize thyself with Reaktoro options and usage, especially when selecting Reaktoro provided databases, database files, and activity models for aqueous, solid, and gas phases. This package does not automatically select any of these options, and will use ideal models by default.**
 
-reaktoro-pse depends on the following packages and/or versions:
+* Please refer here for information on [Reaktoro supported databases](https://reaktoro.org/tutorials/basics/loading-databases.html) (all are supported) 
+* Please refer here for information on [Reaktoro supported activity models](https://reaktoro.org/tutorials/basics/specifying-activity-models.html) (all are supported, included chain operations, or passing in pre-configured activity models)
+
+**By default the package will use:**
+
+* **Database:** [PhreeqcDatabase](https://reaktoro.org/api/classReaktoro_1_1PhreeqcDatabase.html) 
+* **Data file:** [pitzer.dat](https://reaktoro.org/api/classReaktoro_1_1PhreeqcDatabase.html) 
+* **Aqueous activity models:** [ActivityModelIdealAqueous](https://reaktoro.org/api/namespaceReaktoro.html#ae431d4c8a1f283910ae1cf35024091b8)
+* **Gas activity models:** [ActivityModelIdealGas](https://reaktoro.org/api/namespaceReaktoro.html#a7a0788a5a863d987a88b81303d80b427)
+* **Solid activity models:** [ActivityModelIdealSolution](https://reaktoro.org/api/namespaceReaktoro.html#a6581d5c0cde36cae6d9c46dbc32d56f8)
+* **Ion exchange activity modells:** [ActivityModelIonExchange](https://reaktoro.org/api/namespaceReaktoro.html#a6581d5c0cde36cae6d9c46dbc32d56f8)
+
+## 3. Inputs and outputs of the Reaktoro blocks
+The Reaktoro blocks built by this package are designed to solve an equilibrium problem using user provided apparent species or true species, temperature, pressure, and pH, which are broken down to base elements and equilibrated within Reaktoro to provide exact speciation and equilibrium state. Using this state the block can return various information supported by Reaktoro:
+
+* [Chemical properties](https://reaktoro.org/api/classReaktoro_1_1ChemicalProps.html)
+* [Aqueous properties](https://reaktoro.org/api/classReaktoro_1_1AqueousProps.html)
+* Pyomo build properties, which are custom properties built in Pyomo that use chemical properties or aqueous properties as inputs 
+
+Note: Only Reaktoro properties that return a single floating point or real value are supported, as they will be directly matched to a [Pyomo Var](https://pyomo.readthedocs.io/en/stable/pyomo_modeling_components/Variables.html). Reaktoro functions that return arrays or strings are not supported.
+
+## 4. Tutorials, Examples, and Comparisons
+Currently, repo includes several tutorials and examples.
+
+*Tutorials:*
+
+1. [Demonstration of working with Reaktoro block that shows](https://github.com/avdudchenko/reaktoro-pse/blob/main/src/reaktoro_pse/tutorials/ReaktoroBlock_tutorial.ipynb)
+    * How to balance feed charge with ReaktoroBlock
+    * How to add indexed ReaktoroBlocks [WaterTAP RO1D model](https://watertap.readthedocs.io/en/stable/technical_reference/unit_models/reverse_osmosis_1D.html) for calculation of Osmotic pressure
+
+*Examples:*
+
+1. [Example of adding ReaktoroBlock to basic desalination problem](https://github.com/avdudchenko/reaktoro-pse/blob/main/src/reaktoro_pse/examples/simple_desalination.py) that demonstrates how to:
+    * Setup up basic ReaktoroBlock
+    * Calculate basic properties for Scaling Tendency, pH, and Osmotic pressure 
+    * Optimize system pH for operation at target Scaling Tendency
+2. [Example of thermal precipitation](https://github.com/avdudchenko/reaktoro-pse/blob/main/src/reaktoro_pse/examples/thermal_precipitation.py) that demonstrates how to:
+    * Configure different database from default
+    * Get enthalpy and vapor pressure from Reaktoro
+    * Setup precipitation calculation 
+    * Setup simulation for removal of Calcite over different temperatures and estimate required energy input
+3. [Example of ion exchange calculations](https://github.com/avdudchenko/reaktoro-pse/blob/main/src/reaktoro_pse/examples/simple_ion_exchange.py) that demonstrates how to:
+    * Set up ReaktoroBlock for charge neutralizing the feed composition 
+    * Use outputs from speciation block as inputs into a second property block
+    * Add Ion Exchange phase and species into ReaktoroBlock
+    * Optimize addition of acid and bases for maximizing Calcium removal selectivity over Magnesium
+4. [Example of biogas combustion](https://github.com/avdudchenko/reaktoro-pse/blob/main/src/reaktoro_pse/examples/biogas_combustion.py) that demonstrates how to:
+    * Set up ReaktoroBlock for customized database
+    * Use Condensed phase
+
+*Comparisons of Reaktoro-pse to PhreeqC when simulating*
+
+These comparisons further demonstrate how to setup Reaktoro-pse for each type of calculation.
+
+1. [Water removal from solution (e.g evaporative processes)](https://github.com/avdudchenko/reaktoro-pse/blob/main/src/reaktoro_pse/examples/reaktoro_pse_to_phreeqc_comparison/water_removal_comparison.py)
+2. [Vapor pressure calculation at different temperatures](https://github.com/avdudchenko/reaktoro-pse/blob/main/src/reaktoro_pse/examples/reaktoro_pse_to_phreeqc_comparison/vapor_pressure_comparison.py)
+3. [Precipitation of mineral phases](https://github.com/avdudchenko/reaktoro-pse/blob/main/src/reaktoro_pse/examples/reaktoro_pse_to_phreeqc_comparison/precipitation_comparison.py)
+4.  [Mixing of two solution](https://github.com/avdudchenko/reaktoro-pse/blob/main/src/reaktoro_pse/examples/reaktoro_pse_to_phreeqc_comparison/solution_mixing_comparison.py)
+
+## 5. Solvers 
+To-date, this has only been tested with [cyipopt](https://cyipopt.readthedocs.io/en/stable/) - other solvers have not been tested. 
+* For accessing other [HSL linear solvers](http://www.hsl.rl.ac.uk/ipopt/) beyond MUMPS (such as MA27, MA57, etc. solver common to WaterTAP and IDAES-PSE) follow [these instructions](https://cyipopt.readthedocs.io/en/latest/install.html) for adding it to cyipopt.
+
+## 6. Known issues
+
+### Closing dual infeasibility 
+In some cases Ipopt might struggle to close unscaled dual infeasibility even when primal is reduced to <1e-8. This will be typically seen in the solver trace, that shows **inf_pr** being reduced to <1e-8 and **inf_du** stops being reduced and solver does not terminate or terminates with **Solved to Acceptable level**. [This is a result of using approximated hessian in the GrayBox model causing issues in calculations of dual infeasibility error.](https://list.coin-or.org/pipermail/ipopt/2007-February/000700.html)
+
+### Solutions
+A. Set Ipopt solver option "recalc_y" to "yes"
+
+This option will force Ipopt to use least squares method to calculate dual infeasibility, potentially improving accuracy of its estimates and reduce it to below tolerance level. Details on the [recalc_y](https://coin-or.github.io/Ipopt/OPTIONS.html#OPT_recalc_y) and [recalc_y_feas_tol](https://coin-or.github.io/Ipopt/OPTIONS.html#OPT_recalc_y_feas_tol). 
+
+    cy_solver = get_solver(solver="cyipopt-watertap")
+    cy_solver.options['recalc_y']='yes'
+
+
+B. Use exact derivatives instead of numeric
+
+The numeric derivatives carry additional errors that reduce accuracy in estimates of dual infeasibility. You can check which outputs in your Reaktoro block are exact or numeric by using **your_reaktor_block.display_jacobian_outputs()**. 
+
+If option "A" did not work, using exact derivatives can potentially solve this issue. This can be accomplished by using properties with exact derivatives listed in [JacoibanRows class](https://github.com/avdudchenko/reaktoro-pse/blob/868efe883dbc26654b53a32e5a58e8b6ee2af5c7/src/reaktoro_pse/core/reaktoro_jacobian.py#L51). These properties can be used to write Pyomo constraints that calculate the desired property. Some properties are already supported and examples are shown of how to build them in [PyomoProperties](https://github.com/avdudchenko/reaktoro-pse/blob/868efe883dbc26654b53a32e5a58e8b6ee2af5c7/src/reaktoro_pse/core/reaktoro_outputs.py#L118) class. 
+
+Supported PyomoProperties with exact derivatives:
+
+- scalingTendencyDirect - this only designed to work with PhreeqC data bases 
+- phDirect
+- osmoticPressure
+- vaporPressure
+
+These properties are accessed as any other property in ReaktoroBlock. Simply pass ('scalingTendencyDirect',phase) to outputs.  
+
+
+## 7. Requesting new features or issues
+Please include a minimal example using Reaktoro for your specific feature or issue request if possible. Please also include full traceback for errors the occur. 
+
+## 8. Compatibility
+Reaktoro-pse depends on the following packages and/or versions:
 
 - Python 3.9 through 3.12
-- Reaktoro 2.12
+- Reaktoro>=2.12.3
 - CyIpopt 1.4.1
-- Pyomo
+- Pyomo>=6.8.0
+- idaes-pse>=2.5.0
+- watertap>=1.0.0 - (required for watertap-cyipopt wrapper only)
 
-## Getting started (for contributors)
+## 9. Getting started (for contributors)
 
 ### Prerequisites
 
@@ -21,7 +123,7 @@ reaktoro-pse depends on the following packages and/or versions:
 ```sh
 git clone https://github.com/watertap-org/reaktoro-pse.git
 cd reaktoro-pse
-conda create --yes -c conda-forge --name reaktoro-pse-dev python=3.11 reaktoro=2.12.1 cyipopt=1.4.1
+conda create --yes -c conda-forge --name reaktoro-pse-dev python=3.11 reaktoro=2.12.3 cyipopt=1.4.1
 conda activate reaktoro-pse-dev
 pip install -r requirements-dev.txt
 ```

pyproject.toml

@@ -15,6 +15,7 @@ dependencies = [
 [project.optional-dependencies]
 testing = [
     "pytest>=8",
+    "nbmake",
 ]
 
 [tool.setuptools_scm]

requirements-dev.txt

@@ -1,3 +1,6 @@
 black==24.8.0
 pytest-cov
+idaes-pse>=2.5.0
+watertap>=1.0.0
+pyomo>=6.8.0
 --editable .[testing]

src/reaktoro_pse/core/pyomo_property_writer/property_functions.py

@@ -0,0 +1,140 @@
+#################################################################################
+# WaterTAP Copyright (c) 2020-2024, The Regents of the University of California,
+# through Lawrence Berkeley National Laboratory, Oak Ridge National Laboratory,
+# National Renewable Energy Laboratory, and National Energy Technology
+# Laboratory (subject to receipt of any required approvals from the U.S. Dept.
+# of Energy). All rights reserved.
+#
+# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and license
+# information, respectively. These files are also available online at the URL
+# "https://github.com/watertap-org/reaktoro-pse/"
+#################################################################################
+from pyomo.environ import log10, log, exp
+
+
+def build_scaling_tendency_constraint(rkt_output_object):
+    user_output_var = rkt_output_object.pyomo_var
+    build_properties = rkt_output_object.pyomo_build_options.properties
+    return (
+        user_output_var
+        == 10
+        ** build_properties[
+            ("saturationIndex", rkt_output_object.property_index)
+        ].pyomo_var
+    )
+
+
+def build_ph_constraint(rkt_output_object):
+    user_output_var = rkt_output_object.pyomo_var
+    build_properties = rkt_output_object.pyomo_build_options.properties
+    return (
+        -build_properties[("speciesActivityLn", "H+")].pyomo_var / log(10)
+        == user_output_var
+    )
+
+
+def build_vapor_pressure_constraint(rkt_output_object):
+    user_output_var = rkt_output_object.pyomo_var
+    build_properties = rkt_output_object.pyomo_build_options.properties
+    return (
+        exp(
+            build_properties[
+                ("speciesActivityLn", rkt_output_object.property_index)
+            ].pyomo_var
+        )
+        * 101325
+        == user_output_var
+    )
+
+
+# work around provided in https://github.com/reaktoro/reaktoro/discussions/398
+# can not be implemented as reference chemical potential is function of p and t
+# def build_vapor_pressure_direct_constraint(rkt_output_object):
+#     user_output_var = rkt_output_object.pyomo_var
+#     build_properties = rkt_output_object.pyomo_build_options.properties
+#     build_options = rkt_output_object.pyomo_build_options.options
+#     return (
+#         exp(
+#             build_properties[
+#                 ("speciesChemicalPotential", rkt_output_object.property_index)
+#             ].pyomo_var
+#             - build_options["reference_chemical_potential"]
+#         )
+#         / (
+#             build_options["gas_constant"]
+#             * build_properties[("temperature", None)].pyomo_var
+#         )
+#         * 1e5
+#         == user_output_var
+#     )
+
+
+def build_osmotic_constraint(rkt_output_object):
+    # reference  https://help.syscad.net/PHREEQC_Reverse_Osmosis
+    user_output_var = rkt_output_object.pyomo_var
+    build_properties = rkt_output_object.pyomo_build_options.properties
+    build_options = rkt_output_object.pyomo_build_options.options
+    return (
+        user_output_var
+        == -(
+            build_options["gas_constant"]
+            * build_properties[("temperature", None)].pyomo_var
+        )
+        / build_properties[
+            ("speciesStandardVolume", rkt_output_object.property_index)
+        ].pyomo_var
+        * build_properties[
+            ("speciesActivityLn", rkt_output_object.property_index)
+        ].pyomo_var
+    )
+
+
+def build_direct_scaling_tendency_constraint(rkt_output_object):
+    # https://reaktoro.org/api/namespaceReaktoro.html#a55b9a29cdf35e98a6b07e67ed2edbc25
+    user_output_var = rkt_output_object.pyomo_var
+    build_properties = rkt_output_object.pyomo_build_options.properties
+    build_options = rkt_output_object.pyomo_build_options.options
+    # TODO: Needs to add temperature unit verificaiton and pressure unit verification
+    temperature_var = build_properties[("temperature", None)].pyomo_var
+    if build_options["logk_type"] == "Analytical":
+        A_params = build_options["logk_paramters"]
+        log_k = [A_params["A1"]]
+        # temp dependence for phreeqc
+        if A_params["A2"] != 0:
+            log_k.append(A_params["A2"] * temperature_var)
+        if A_params["A3"] != 0:
+            log_k.append(A_params["A3"] * temperature_var**-1)
+        if A_params["A4"] != 0:
+            log_k.append(A_params["A4"] * log10(temperature_var))
+        if A_params["A5"] != 0:
+            log_k.append(A_params["A5"] * temperature_var**-2)
+        if A_params["A6"] != 0:
+            log_k.append(A_params["A6"] * temperature_var**2)
+
+    if build_options["logk_type"] == "VantHoff":
+        vfparams = build_options["logk_paramters"]
+        log_k = [
+            vfparams["lgKr"]
+            - vfparams["dHr"]
+            / build_options["gas_constant"]
+            * (1 / temperature_var - 1 / vfparams["Tr"])
+        ]
+    # pressure dependenance
+    log_k.append(
+        -(
+            build_options["delta_V"]
+            * (build_properties[("pressure", None)].pyomo_var - 101325)
+            / (log(10) * build_options["gas_constant"] * temperature_var)
+        )
+    )
+
+    activities = []
+    for key, obj in build_properties.items():
+        if "speciesActivityLn" in key:
+            activities.append(obj.pyomo_var * obj.stoichiometric_coeff / log(10))
+    return user_output_var == 10 ** (
+        sum(activities)
+        + sum(
+            log_k
+        )  # this is postive here and in log10 fom, so we add instead of subtract
+    )