Performance Improvement of erfint

CHANGELOG.md

@@ -17,6 +17,10 @@
 * [Issue 229](https://github.com/MassimoCimmino/pygfunction/issues/229), [Issue 247](https://github.com/MassimoCimmino/pygfunction/issues/247) - Added citation to IGSHPA conference paper on *pygfunction* v2.2 in the documention. Added a `CITATION.cff` file to suggest a correct citation on github.
 * [Issue 230](https://github.com/MassimoCimmino/pygfunction/issues/230) - Configured github actions to publish *pygfunction* on Pypi on creation of a release on github.
 
+### Enhancements
+
+* [Issue 204](https://github.com/MassimoCimmino/pygfunction/issues/204) - Added support for Python 3.9 and 3.10. [CoolProp](https://www.coolprop.org/) is removed from the dependencies and replace with [SecondaryCoolantProps](https://github.com/mitchute/SecondaryCoolantProps).
+
 ## Version 2.2.1 (2022-08-12)
 
 ### Bug fixes

pygfunction/heat_transfer.py

@@ -4,7 +4,7 @@
 from scipy.special import erfc, erf, roots_legendre
 
 from .boreholes import Borehole
-from .utilities import erfint, exp1, _erf_coeffs
+from .utilities import erf_int, exp1, _erf_coeffs
 
 
 def finite_line_source(
@@ -1119,7 +1119,7 @@ def _finite_line_source_integrand(dis, H1, D1, H2, D2, reaSource, imgSource):
                       D2 + D1 + H1,
                       D2 + D1 + H2 + H1],
                      axis=-1)
-        f = lambda s: s**-2 * np.exp(-dis**2*s**2) * np.inner(p, erfint(q*s))
+        f = lambda s: s**-2 * np.exp(-dis**2*s**2) * np.inner(p, erf_int(q * s))
     elif reaSource:
         # Real FLS solution
         p = np.array([1, -1, 1, -1])
@@ -1128,7 +1128,7 @@ def _finite_line_source_integrand(dis, H1, D1, H2, D2, reaSource, imgSource):
                       D2 - D1 - H1,
                       D2 - D1 + H2 - H1],
                      axis=-1)
-        f = lambda s: s**-2 * np.exp(-dis**2*s**2) * np.inner(p, erfint(q*s))
+        f = lambda s: s**-2 * np.exp(-dis**2*s**2) * np.inner(p, erf_int(q * s))
     elif imgSource:
         # Image FLS solution
         p = np.array([1, -1, 1, -1])
@@ -1137,7 +1137,7 @@ def _finite_line_source_integrand(dis, H1, D1, H2, D2, reaSource, imgSource):
                       D2 + D1 + H1,
                       D2 + D1 + H2 + H1],
                      axis=-1)
-        f = lambda s: s**-2 * np.exp(-dis**2*s**2) * np.inner(p, erfint(q*s))
+        f = lambda s: s**-2 * np.exp(-dis**2*s**2) * np.inner(p, erf_int(q * s))
     else:
         # No heat source
         f = lambda s: np.zeros(np.broadcast_shapes(
@@ -1304,7 +1304,7 @@ def _finite_line_source_equivalent_boreholes_integrand(dis, wDis, H1, D1, H2, D2
                       D2 + D1 + H1,
                       D2 + D1 + H2 + H1],
                      axis=-1)
-        f = lambda s: s**-2 * (np.exp(-dis**2*s**2) @ wDis).T * np.inner(p, erfint(q*s))
+        f = lambda s: s**-2 * (np.exp(-dis**2*s**2) @ wDis).T * np.inner(p, erf_int(q * s))
     elif reaSource:
         # Real FLS solution
         p = np.array([1, -1, 1, -1])
@@ -1313,7 +1313,7 @@ def _finite_line_source_equivalent_boreholes_integrand(dis, wDis, H1, D1, H2, D2
                       D2 - D1 - H1,
                       D2 - D1 + H2 - H1],
                      axis=-1)
-        f = lambda s: s**-2 * (np.exp(-dis**2*s**2) @ wDis).T * np.inner(p, erfint(q*s))
+        f = lambda s: s**-2 * (np.exp(-dis**2*s**2) @ wDis).T * np.inner(p, erf_int(q * s))
     elif imgSource:
         # Image FLS solution
         p = np.array([1, -1, 1, -1])
@@ -1322,7 +1322,7 @@ def _finite_line_source_equivalent_boreholes_integrand(dis, wDis, H1, D1, H2, D2
                       D2 + D1 + H1,
                       D2 + D1 + H2 + H1],
                      axis=-1)
-        f = lambda s: s**-2 * (np.exp(-dis**2*s**2) @ wDis).T * np.inner(p, erfint(q*s))
+        f = lambda s: s**-2 * (np.exp(-dis**2*s**2) @ wDis).T * np.inner(p, erf_int(q * s))
     else:
         # No heat source
         f = lambda s: np.zeros(np.broadcast_shapes(

pygfunction/media.py

@@ -192,4 +192,4 @@ def Prandlt_number(self):
             Prandlt number.
 
         """
-        return self.fluid.prandtl(self.T_C)
+        return self.fluid.prandtl(self.T_C)

pygfunction/utilities.py

@@ -4,6 +4,8 @@
 import numpy.polynomial.polynomial as poly
 from scipy.special import erf
 import warnings
+from typing import Union
+from numpy.typing import NDArray
 
 
 def cardinal_point(direction):
@@ -22,7 +24,33 @@ def cardinal_point(direction):
     return compass[direction]
 
 
-def erfint(x):
+sqrt_pi = 1 / np.sqrt(np.pi)
+
+
+def erf_int(x: Union[NDArray[np.float64], float]) -> NDArray[np.float64]:
+    """
+    Integral of the error function.
+
+    Parameters
+    ----------
+    x : float or array
+        Argument.
+
+    Returns
+    -------
+    float or array
+        Integral of the error function.
+
+    """
+    abs_x = np.abs(x)
+    y_new = abs_x-sqrt_pi
+    idx = np.less(abs_x, 4)
+    abs_2 = abs_x[idx]
+    y_new[idx] = abs_2 * erf(abs_2) - (1.0 - np.exp(-abs_2*abs_2)) * sqrt_pi
+    return y_new
+
+
+def erf_int_old(x: Union[NDArray[np.float64], float]) -> NDArray[np.float64]:
     """
     Integral of the error function.
 
@@ -37,7 +65,7 @@ def erfint(x):
         Integral of the error function.
 
     """
-    return x * erf(x) - 1.0 / np.sqrt(np.pi) * (1.0 - np.exp(-x**2))
+    return x * erf(x) - (1.0 - np.exp(-x*x)) / np.sqrt(np.pi)
 
 
 def exp1(x):

pyproject.toml

@@ -3,7 +3,7 @@ requires = ["setuptools", "wheel"]
 build-backend = "setuptools.build_meta"
 
 [tool.pytest.ini_options]
-addopts = "--cov=pygfunction"
+# addopts = "--cov=pygfunction"
 testpaths = [
     "tests",
 ]

tests/test_erf_int.py

@@ -0,0 +1,45 @@
+from pygfunction.utilities import erf_int, erf_int_old
+import numpy as np
+from time import perf_counter_ns
+
+
+def test_erf_int_error():
+    print(f'test of erf error and performance')
+    x = np.arange(-3.9, 3.9, 0.01)
+    tic = perf_counter_ns()
+    y_new = erf_int(x)
+    toc = perf_counter_ns()
+    dt_new1 = toc-tic
+    tic = perf_counter_ns()
+    y = erf_int_old(x)
+    toc = perf_counter_ns()
+    dt_old1 = toc - tic
+    assert np.allclose(y, y_new, rtol=1e-10)
+
+    print(f'new time {dt_new1 / 1_000_000} ms; old time { dt_old1 / 1_000_000} ms')
+
+    x = np.arange(-500, 500, 5)
+    tic = perf_counter_ns()
+    y_new = erf_int(x)
+    toc = perf_counter_ns()
+    dt_new2 = toc - tic
+    tic = perf_counter_ns()
+    y = erf_int_old(x)
+    toc = perf_counter_ns()
+    dt_old2 = toc - tic
+    assert np.allclose(y, y_new, rtol=1e-10)
+
+    print(f'new time {dt_new2 / 1_000_000} ms; old time {dt_old2 / 1_000_000} ms')
+
+    x = np.arange(-500, 500, 0.01)
+    tic = perf_counter_ns()
+    y_new = erf_int(x)
+    toc = perf_counter_ns()
+    dt_new3 = toc - tic
+    tic = perf_counter_ns()
+    y = erf_int_old(x)
+    toc = perf_counter_ns()
+    dt_old3 = toc - tic
+    assert np.allclose(y, y_new, rtol=1e-10)
+
+    print(f'new time {dt_new3/1_000_000} ms; old time {dt_old3/1_000_000} ms')