adding library for nonlinear watermodels

Physics/water/LICENSE

+Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification,
+are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+2. Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+3. Neither the name of the OpenGeoSys Community nor the
+   names of its contributors may be used to endorse or promote products
+   derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
+INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
+NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
+OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
+OF SUCH DAMAGE.
+
+

Physics/water/README.md

+# water.py
+class for creating nonlinear watermodels to be used in OGS
+Most data taken from https://www.engineeringtoolbox.com
\ No newline at end of file

Physics/water/setup.py

+# -*- coding: utf-8 -*-
+from setuptools import setup
+
+setup(name="water",
+      version=0.01,
+      maintainer="Jörg Buchwald",
+      maintainer_email="joerg_buchwald@ufz.de",
+      author="Jörg Buchwald",
+      author_email="joerg.buchwald@ufz.de",
+      url="https://github.com/joergbuchwald/mini-projects",
+      classifiers=["Intended Audience :: Science/Research",
+          "Topic :: Scientific/Engineering :: Visualization",
+          "Topic :: Scientific/Engineering :: Physics",
+          "Topic :: Scientific/Engineering :: Mathematics",
+          "License :: OSI Approved :: MIT License",
+          "Programming Language :: Python :: 3",
+          "Programming Language :: Python :: 3.8"],
+      license="BSD-3 -  see LICENSE.txt",
+      platforms=["Windows", "Linux", "Solaris", "Mac OS-X", "Unix"],
+      include_package_data=True,
+      python_requires='>=3.8',
+      install_requires=["numpy"],
+      py_modules=["water/water"],
+      packages=["water/properties"])

Physics/water/water/properties/conductivity.py

+"""
+Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+              Distributed under a Modified BSD License.
+                See accompanying file LICENSE or
+                http://www.opengeosys.org/project/license
+
+"""
+# pylint: disable=C0103, R0902, R0914, R0913
+import numpy as np
+from water.properties import template
+T0 = 273.15
+
+class conductivity_1(template.PROPERTY):
+    def value(self, temperature):
+        # source: https://www.engineeringtoolbox.com/water-liquid-gas-thermal-conductivity-temperature-pressure-d_2012.html
+        temp = np.array([0.01, 10., 20., 30., 40., 50., 60., 70., 80., 90., 99.6])+T0
+        if (np.min(temperature) < temp[0]) or (np.max(temperature) > temp[-1]):
+            print("The temperature is not within the defined domain")
+        K = np.array([
+            0.556, 0.579, 0.598, 0.614, 0.629, 0.641, 0.651, 0.660, 0.667,
+            0.673, 0.677
+        ])
+        return np.interp(temperature, temp, K)
+    def dvalue(self, temperature):
+        return np.gradient(self.value(temperature),temperature)
+    def exprtk_value(self):
+        raise NotImplementedError
+    def exprtk_dvalue(self):
+        raise NotImplementedError
+
+class conductivity_2(template.PROPERTY):
+    def value(self, temperature):
+        # use only one-liner to keep it parsable:
+        return  3.65759470e-08 * temperature**3 - 4.51285141e-05 * temperature**2 + 1.88275537e-02 * temperature - 1.96475421e+00
+    def dvalue(self, temperature):
+        # use only one-liner to keep it parsable:
+        return 3.65759470e-08 * 3 * temperature**2 - 4.51285141e-05 * 2 * temperature + 1.88275537e-02
+    def dvaluenum(self,temperature):
+        return np.gradient(self.value(temperature),temperature)
+    def exprtk_value(self):
+        string = self._getcodeasstring("value", None)
+        string = self._convertpythontoexprtk(string)
+        return string
+    def exprtk_dvalue(self):
+        string = self._getcodeasstring("dvalue",None)
+        string = self._convertpythontoexprtk(string)
+        return string

Physics/water/water/properties/density.py

+"""
+Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+              Distributed under a Modified BSD License.
+                See accompanying file LICENSE or
+                http://www.opengeosys.org/project/license
+
+"""
+# pylint: disable=C0103, R0902, R0914, R0913
+import numpy as np
+from water.properties import template
+T0 = 273.15
+
+class density_1(template.PROPERTY):
+    def value(self, temperature):
+        # source: https://www.engineeringtoolbox.com/water-density-specific-weight-d_595.html
+        temp = np.array([
+            0.1, 1.0, 4., 10., 15., 20., 25., 30., 35., 40., 45., 50., 60.,
+            70., 80., 90., 100., 110., 140., 200., 260., 300., 360.
+        ]) + T0
+        if (np.min(temperature) < temp[0]) or (np.max(temperature) > temp[-1]):
+            print("The temperature is not within the defined domain")
+        dens = np.array([
+            999.85, 999.9, 999.97, 999.7, 999.1, 998.2, 997., 995.6, 994.,
+            992.2, 990.2, 988., 983.2, 977.76, 971.8, 965.3, 958.6, 951.0,
+            926.1, 865., 783.6, 712.2, 527.6
+        ])
+        return np.interp(temperature, temp, dens)
+    def dvalue(self,temperature):
+        return np.gradient(self.value(temperature),temperature)
+    def exprtk_value(self):
+        raise NotImplementedError
+    def exprtk_dvalue(self):
+        raise NotImplementedError
+
+class density_2(template.PROPERTY):
+    def value(self, temperature, phase_pressure):
+        # use only one-liner to keep it parsable:
+        return 1000.1 * (1 - (-6*10**-6 * (temperature-T0)**4 + 0.001667 * (temperature-T0)**3 + -0.197796 * (temperature-T0)**2 + 16.86446 * (temperature-T0) - 64.319951)/10**6 * (temperature-293.15) + 4.65e-10 * (np.maximum(phase_pressure, 0.0) - 1e5))
+    def dvalue(self, temperature, phase_pressure, variable="temperature"):
+        if variable == "temperature":
+            # use only one-liner to keep it parsable:
+            return -1000.1 * ((-6*10**-6 * (temperature-T0)**4 + 0.001667 * (temperature-T0)**3 + -0.197796 * (temperature-T0)**2 + 16.86446 * (temperature-T0) - 64.319951)/10**6+(-6*10**-6 * 4 * (temperature-T0)**3 + 0.001667 * 3 * (temperature-T0)**2 + -0.197796 * 2 * (temperature-T0) + 16.86446 )/10**6*(temperature-293.15))
+        elif variable == "phase_pressure":
+            # use only one-liner to keep it parsable:
+            return np.maximum(np.sign(phase_pressure),0)* 1000.1 * 4.65e-10
+    def dvaluenum(self,temperature, phase_pressure, variable="temperature"):
+        if variable == "temperature":
+            return np.gradient(self.value(temperature,phase_pressure),temperature)
+        elif variable == "phase_pressure":
+            return np.gradient(self.value(temperature,phase_pressure),phase_pressure)
+    def exprtk_value(self):
+        string = self._getcodeasstring("value", None)
+        string = self._convertpythontoexprtk(string)
+        return string
+    def exprtk_dvalue(self, variable="temperature"):
+        string = self._getcodeasstring("dvalue",variable)
+        string = self._convertpythontoexprtk(string)
+        return string
+
+class density_3(template.PROPERTY):
+    def value(self, temperature, phase_pressure):
+        # use only one-liner to keep it parsable:
+        return (5.38784232e-06*temperature**3 - 8.61928859e-03*temperature**2 +  3.44638154e+00*temperature + 5.92596450e+02) * np.exp(4.65e-10 * (np.maximum(phase_pressure, 0.0) - 1e5))
+    def dvalue(self, temperature, phase_pressure, variable="temperature"):
+        if variable == "temperature":
+            return (3*5.38784232e-06*temperature**2 - 2*8.61928859e-03*temperature +  3.44638154e+00) * np.exp(4.65e-10 * (np.maximum(phase_pressure, 0.0) - 1e5))
+        elif variable == "phase_pressure":
+            return (5.38784232e-06*temperature**3 - 8.61928859e-03*temperature**2 +  3.44638154e+00*temperature + 5.92596450e+02) * np.maximum(np.sign(phase_pressure),0) * 4.65e-10
+    def dvaluenum(self,temperature, phase_pressure, variable="temperature"):
+        if variable == "temperature":
+            # use only one-liner to keep it parsable:
+            return np.gradient(self.value(temperature,phase_pressure),temperature)
+        elif variable == "phase_pressure":
+            # use only one-liner to keep it parsable:
+            return np.gradient(self.value(temperature,phase_pressure),phase_pressure)
+    def exprtk_value(self):
+        string = self._getcodeasstring("value", None)
+        string = self._convertpythontoexprtk(string)
+        return string
+    def exprtk_dvalue(self, variable="temperature"):
+        string = self._getcodeasstring("dvalue",variable)
+        string = self._convertpythontoexprtk(string)
+        return string

Physics/water/water/properties/expansivity.py

+"""
+Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+              Distributed under a Modified BSD License.
+                See accompanying file LICENSE or
+                http://www.opengeosys.org/project/license
+
+"""
+# pylint: disable=C0103, R0902, R0914, R0913
+import numpy as np
+from water.properties import template
+T0 = 273.15
+
+class expansivity_1(template.PROPERTY):
+    def value(self, temperature):
+        # source: https://www.engineeringtoolbox.com/water-density-specific-weight-d_595.html
+        temp = np.array([
+            0.0, 4., 10., 20., 30., 40., 50., 60., 70., 80., 90., 140., 200.,
+            260.
+        ]) + T0
+        if (np.min(temperature) < temp[0]) or (np.max(temperature) > temp[-1]):
+            print("The temperature is not within the defined domain")
+        beta = np.array([
+            -5.e-5, 0.003e-4, 8.8e-5, 2.07e-4, 3.03e-4, 3.85e-4, 4.57e-4, 5.22e-4,
+            5.82e-4, 6.40e-4, 6.95e-4, 9.75e-4, 1.59e-3, 2.21e-3
+        ])
+        return np.interp(temperature, temp, beta)
+    def dvalue(self, temperature):
+        return np.gradient(self.value(temperature), temperature)
+    def exprtk_value(self):
+        raise NotImplementedError
+    def exprtk_dvalue(self):
+        raise NotImplementedError
+
+class expansivity_2(template.PROPERTY):
+    def value(self, temperature):
+        # use only one-liner to keep it parsable:
+        return  -1.46673616e-12 * temperature**4 + 2.46256925e-09 * temperature**3 - 1.51400613e-06 * temperature**2 + 4.11530308e-04 * temperature - 4.15284133e-02
+    def dvalue(self, temperature):
+        # use only one-liner to keep it parsable:
+        return -1.46673616e-12 * 4 * temperature**3 + 2.46256925e-09 * 3 * temperature**2 - 1.51400613e-06 * 2 * temperature + 4.11530308e-04
+    def dvaluenum(self,temperature):
+        return np.gradient(self.value(temperature),temperature)
+    def exprtk_value(self):
+        string = self._getcodeasstring("value", None)
+        string = self._convertpythontoexprtk(string)
+        return string
+    def exprtk_dvalue(self):
+        string = self._getcodeasstring("dvalue",None)
+        string = self._convertpythontoexprtk(string)
+        return string

Physics/water/water/properties/specificheatcapacity.py

+"""
+Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+              Distributed under a Modified BSD License.
+                See accompanying file LICENSE or
+                http://www.opengeosys.org/project/license
+
+"""
+# pylint: disable=C0103, R0902, R0914, R0913
+import numpy as np
+from water.properties import template
+T0 = 273.15
+
+class specificheatcapacity_1(template.PROPERTY):
+    def value(self, temperature):
+        # source: https://www.engineeringtoolbox.com/specific-heat-capacity-water-d_660.html
+        temp = np.array([0.01, 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.]) + T0
+        if (np.min(temperature) < temp[0]) or (np.max(temperature) > temp[-1]):
+            print("The temperature is not within the defined domain")
+        cw = np.array([
+            4217., 4191., 4157., 4118., 4074., 4026., 3977., 3925., 3873.,
+            3820., 3768
+        ])
+        return np.interp(temperature, temp, cw)
+    def dvalue(self, temperature):
+        return np.gradient(self.value(temperature),temperature)
+    def exprtk_value(self):
+        raise NotImplementedError
+    def exprtk_dvalue(self):
+        raise NotImplementedError
+
+class specificheatcapacity_2(template.PROPERTY):
+    def value(self, temperature):
+        # use only one-liner to keep it parsable:
+        return  1.55452794e-04 * temperature**3 - 1.64289282e-01 * temperature**2 + 5.25998084e+01 * temperature - 1.06049924e+03
+    def dvalue(self, temperature):
+        # use only one-liner to keep it parsable:
+        return 1.55452794e-04 * 3 * temperature**2 - 1.64289282e-01 * 2 * temperature + 5.25998084e+01
+    def dvaluenum(self,temperature):
+        return np.gradient(self.value(temperature),temperature)
+    def exprtk_value(self):
+        string = self._getcodeasstring("value", None)
+        string = self._convertpythontoexprtk(string)
+        return string
+    def exprtk_dvalue(self):
+        string = self._getcodeasstring("dvalue",None)
+        string = self._convertpythontoexprtk(string)
+        return string

Physics/water/water/properties/template.py

+import pathlib
+class PROPERTY:
+    def __init__(self):
+        pass
+    def value(self):
+        raise NotImplementedError
+    def dvalue(self):
+        raise NotImplementedError
+    def exprtk_value(self):
+        raise NotImplementedError
+    def exprtk_dvalue(self):
+        raise NotImplementedError
+    def _getcodeasstring(self, fct, variable):
+        classname = self.__class__.__name__
+        string = ""
+        currentpath = pathlib.Path(__file__).parent.resolve()
+        with open(f"{currentpath}/{classname.split('_')[0]}.py") as source:
+            currentclass = None
+            currentmethod = None
+            currentvariable = None
+            for line in source:
+                linelist = line.split(" ")
+                try:
+                    if linelist[0] == "class":
+                        currentclass = linelist[1].split("(")[0]
+                        currentmethod = None
+                        currentvariable = None
+                except IndexError:
+                    pass
+                try:
+                    if linelist[4] == "def":
+                        currentmethod = linelist[5].split("(")[0]
+                except IndexError:
+                    pass
+                try:
+                    if (linelist[8] == "if") or (linelist[8] == "elif"):
+                        if linelist[9] == "variable":
+                            currentvariable = linelist[11].split('"')[1]
+                except IndexError:
+                    pass
+                try:
+                    if classname == currentclass:
+                        if currentmethod == fct:
+                            if currentvariable == variable:
+                                if "return" in linelist:
+                                    string = line.split("return")[1]
+                except IndexError:
+                    pass
+        return string
+
+    def _convertpythontoexprtk(self, string):
+        string = string.replace("**", "^")
+        string = string.replace("np.maximum","max")
+        string = string.replace("np.minimum","min")
+        string = string.replace("np.sign","sgn")
+        string = string.replace("np.","")
+        string = string.replace("\n","")
+        return string

Physics/water/water/properties/viscosity.py

+"""
+Copyright (c) 2012-2021, OpenGeoSys Comunity (http://www.opengeosys.org)
+              Distributed under a Modified BSD License.
+                See accompanying file LICENSE or
+                http://www.opengeosys.org/project/license
+
+"""
+# pylint: disable=C0103, R0902, R0914, R0913
+import numpy as np
+from water.properties import template
+T0 = 273.15
+
+class viscosity_1(template.PROPERTY):
+    def value(self, temperature):
+        # source: https://www.engineeringtoolbox.com/water-dynamic-kinematic-viscosity-d_596.html
+        temp = np.array([10., 20., 30., 40., 50., 60., 70., 80., 90., 100.]) + T0
+        if (np.min(temperature) < temp[0]) or (np.max(temperature) > temp[-1]):
+            print("The temperature is not within the defined domain")
+        vis = np.array([
+            0.0013, 0.001, 0.0007978, 0.0006531, 0.0005471, 0.0004658,
+            0.0004044, 0.000355, 0.000315, 0.00002822
+        ])
+        return np.interp(temperature, temp, vis)
+    def dvalue(self, temperature):
+        return np.gradient(self.value(temperature),temperature)
+    def exprtk_value(self):
+        raise NotImplementedError
+    def exprtk_dvalue(self):
+        raise NotImplementedError
+
+class viscosity_2(template.PROPERTY):
+    def value(self, temperature):
+        # use only one-liner to keep it parsable:
+        return -2.01570959e-09 * temperature**3 + 2.11402803e-06 * temperature**2 - 7.43932150e-04 * temperature +  8.82066680e-02
+    def dvalue(self, temperature):
+        # use only one-liner to keep it parsable:
+        return -2.01570959e-09 * 3 * temperature**2 + 2.11402803e-06 * 2 * temperature - 7.43932150e-04
+    def dvaluenum(self,temperature):
+        return np.gradient(self.value(temperature),temperature)
+    def exprtk_value(self):
+        string = self._getcodeasstring("value", None)
+        string = self._convertpythontoexprtk(string)
+        return string
+    def exprtk_dvalue(self):
+        string = self._getcodeasstring("dvalue",None)
+        string = self._convertpythontoexprtk(string)
+        return string
+
+class viscosity_3(template.PROPERTY):
+    def value(self, temperature):
+        return 2.414*10**(-5)*10**(247.8/(temperature-140))
+    def dvalue(self, temperature):
+        return -2.414*10**(-5)*247.8*np.log(10) * 10**(247.8/(temperature-140))/(temperature-140)**2
+    def dvaluenum(self,temperature):
+        return np.gradient(self.value(temperature),temperature)
+    def exprtk_value(self):
+        string = self._getcodeasstring("value", None)
+        string = self._convertpythontoexprtk(string)
+        return string
+    def exprtk_dvalue(self):
+        string = self._getcodeasstring("dvalue",None)
+        string = self._convertpythontoexprtk(string)
+        return string
+
+class viscosity_4(template.PROPERTY):
+    def value(self, temperature):
+        return 4.2844e-5+1/(0.157 * (temperature-208.157) * (temperature-208.157)-91.296)
+    def dvalue(self, temperature):
+        raise NotImplementedError
+    def dvaluenum(self,temperature):
+        return np.gradient(self.value(temperature),temperature)
+    def exprtk_value(self):
+        string = self._getcodeasstring("value", None)
+        string = self._convertpythontoexprtk(string)
+        return string
+    def exprtk_dvalue(self):
+        string = self._getcodeasstring("dvalue", None)
+        string = self._convertpythontoexprtk(string)
+        return string

Physics/water/water/water.py

+"""
+Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+              Distributed under a Modified BSD License.
+                See accompanying file LICENSE or
+                http://www.opengeosys.org/project/license
+
+"""
+
+# pylint: disable=C0103, R0902, R0914, R0913
+from water.properties import (density, viscosity, conductivity, specificheatcapacity, expansivity)
+
+class WATER:
+    """
+    class defining water related properties at atmospheric
+    and arbitrary pressure
+    """
+    def __init__(self, rho=1, mu=0, K=0, c=0, a=0):
+        self.densities = [density.density_1,
+                density.density_2,
+                density.density_3]
+        self.viscosities = [viscosity.viscosity_1,
+                viscosity.viscosity_2,
+                viscosity.viscosity_3,
+                viscosity.viscosity_4]
+        self.conductivities = [conductivity.conductivity_1, conductivity.conductivity_2]
+        self.specificheatcapacities = [specificheatcapacity.specificheatcapacity_1, specificheatcapacity.specificheatcapacity_2]
+        self.expansivities = [expansivity.expansivity_1, expansivity.expansivity_2]
+        self.rho = self.densities[rho]()
+        self.mu = self.viscosities[mu]()
+        self.K = self.conductivities[K]()
+        self.c = self.specificheatcapacities[c]()
+        self.a = self.expansivities[a]()