new functionality tests for scallib (pytest)

python/scallib001/tests/README.md

@@ -0,0 +1,15 @@
+## Collection of functionality and cpu performance tests for scallib001
+
+The tests make use of the standard python test library [pytest](https://docs.pytest.org/en/7.1.x/getting-started.html). 
+
+When necessary, pytest can be installed using pip: pip install pytest.
+
+If in addition you want to do cpu performance benchmarking, you need to
+install the pytest plugin [pytest-benchmark](https://pytest-benchmark.readthedocs.io/en/latest/).
+
+The pytest-benchmark plugin for pytest can be installed using pip: pip install pytest-benchmark.
+
+To run the tests, go to the directory that contains scallib001 and do
+
+python -m pytest scallib001
+

python/scallib001/tests/conftest.py

@@ -0,0 +1,14 @@
+#------------------------------------------------------------------------------------------------
+#  Copyright (c) Shell Global Solutions International B.V. All rights reserved.
+#  Licensed under the MIT License. See License.txt in the project root for license information.
+#------------------------------------------------------------------------------------------------
+
+import pytest
+
+try:
+   import pytest_benchmark
+except:
+   # pytest_benchmark does not installed, skip it
+   @pytest.fixture
+   def benchmark(*args, **kwargs):
+       return None

python/scallib001/tests/test_power_eps1.py

@@ -0,0 +1,208 @@
+#------------------------------------------------------------------------------------------------
+#  Copyright (c) Shell Global Solutions International B.V. All rights reserved.
+#  Licensed under the MIT License. See License.txt in the project root for license information.
+#------------------------------------------------------------------------------------------------
+
+import numpy as np
+import scallib001.relpermlib001 as rlplib
+
+
+def test_power_eps1():
+
+    eps = 1e-4
+
+    # check power_eps(0.0, 0.0, eps), i.e., zero power at zero
+
+    exact_fd_check1 = 1/eps / np.log(1 + 1/eps)
+
+    assert np.allclose( rlplib.power_eps(0.0, 0.0, eps), [0.0, exact_fd_check1] )
+    assert np.allclose( rlplib.power_eps(0.0, 1e-10, eps), [0.0, exact_fd_check1] )
+
+    # check power_eps(x, 0.0, eps), i.e., zero power at zero and finite x
+
+    sv2 = np.logspace(-6, 0, 51)
+
+    exact_f_check2 = np.log(1.0 + sv2/eps)/np.log(1.0 + 1/eps)
+    exact_fd_check2 = 1/(1.0 + sv2/eps)/eps/np.log(1.0 + 1/eps)
+
+    f2, fd2 = rlplib.power_eps(sv2, 0.0, eps)
+
+    assert np.allclose(f2, exact_f_check2)
+    assert np.allclose(fd2, exact_fd_check2)
+
+
+
+def test_power_eps2():
+
+    eps = 1e-6
+
+    sv1 = np.logspace(-6, 0, 51)
+
+    # check that numba compiled version equals python version for N = 0
+
+    numba_f1, numba_fd1 = rlplib.power_eps( sv1, 0.0, eps)
+    python_f1, python_fd1 = rlplib.power_eps.py_func( sv1, 0.0, eps)
+
+    assert np.allclose(numba_f1, python_f1)
+    assert np.allclose(numba_fd1, python_fd1)
+
+    # check that numba compiled version equals python version for N = 0.3
+
+    numba_f2, numba_fd2 = rlplib.power_eps( sv1, 0.3, eps)
+    python_f2, python_fd2 = rlplib.power_eps.py_func( sv1, 0.3, eps)
+
+    assert np.allclose(numba_f2, python_f2)
+    assert np.allclose(numba_fd2, python_fd2)
+
+    # check that numba compiled version equals python version for N = 1.0
+
+    numba_f3, numba_fd3 = rlplib.power_eps( sv1, 1.0, eps)
+    python_f3, python_fd3 = rlplib.power_eps.py_func( sv1, 1.0, eps)
+
+    assert np.allclose(numba_f3, python_f3)
+    assert np.allclose(numba_fd3, python_fd3)
+
+    # check that numba compiled version equals python version for N = 3.1
+
+    numba_f4, numba_fd4 = rlplib.power_eps( sv1, 3.1, eps)
+    python_f4, python_fd4 = rlplib.power_eps.py_func( sv1, 3.1, eps)
+
+    assert np.allclose(numba_f4, python_f4)
+    assert np.allclose(numba_fd4, python_fd4)
+
+def test_power_eps3():
+
+    # Take a very small value for eps
+    eps = 1e-12
+
+    sv1 = np.logspace(-6, 0, 51)
+
+    # check that numba compiled version equals python version for N = 0
+
+    numba_f1, numba_fd1 = rlplib.power_eps( sv1, 0.0, eps)
+    python_f1, python_fd1 = rlplib.power_eps.py_func( sv1, 0.0, eps)
+
+    assert np.allclose(numba_f1, python_f1)
+    assert np.allclose(numba_fd1, python_fd1)
+
+    # check that numba compiled version equals python version for N = 0.3
+
+    numba_f2, numba_fd2 = rlplib.power_eps( sv1, 0.3, eps)
+    python_f2, python_fd2 = rlplib.power_eps.py_func( sv1, 0.3, eps)
+
+    assert np.allclose(numba_f2, python_f2)
+    assert np.allclose(numba_fd2, python_fd2)
+
+    # check that numba compiled version equals python version for N = 1.0
+
+    numba_f3, numba_fd3 = rlplib.power_eps( sv1, 1.0, eps)
+    python_f3, python_fd3 = rlplib.power_eps.py_func( sv1, 1.0, eps)
+
+    assert np.allclose(numba_f3, python_f3)
+    assert np.allclose(numba_fd3, python_fd3)
+
+    # check that numba compiled version equals python version for N = 3.1
+
+    numba_f4, numba_fd4 = rlplib.power_eps( sv1, 3.1, eps)
+    python_f4, python_fd4 = rlplib.power_eps.py_func( sv1, 3.1, eps)
+
+    assert np.allclose(numba_f4, python_f4)
+    assert np.allclose(numba_fd4, python_fd4)
+
+def test_power_eps4():
+
+    eps = 0.0
+
+    sv1 = np.logspace(-6, 0, 51)
+
+    # check that for eps=0 power_eps equals standard np.power for N = 0.0
+
+    numba_f1, numba_fd1 = rlplib.power_eps( sv1, 0.0, eps)
+    numpy_f1, numpy_fd1 = np.power( sv1, 0.0), np.power(sv1, -1.0)*0.0
+
+    assert np.allclose(numba_f1, numpy_f1)
+    assert np.allclose(numba_fd1, numpy_fd1)
+
+    # check that for eps=0 power_eps equals standard np.power for N = 0.3
+
+    numba_f2, numba_fd2 = rlplib.power_eps( sv1, 0.3, eps)
+    numpy_f2, numpy_fd2 = np.power( sv1, 0.3), np.power( sv1, 0.3-1.0 )*0.3
+
+    assert np.allclose(numba_f2, numpy_f2)
+    assert np.allclose(numba_fd2, numpy_fd2)
+
+    # check that for eps=0 power_eps equals standard np.power for N = 1.0
+
+    numba_f3, numba_fd3 = rlplib.power_eps( sv1, 1.0, eps)
+    numpy_f3, numpy_fd3 = np.power( sv1, 1.0), np.power( sv1, 1.0-1.0 )*1.0
+
+    assert np.allclose(numba_f3, numpy_f3)
+    assert np.allclose(numba_fd3, numpy_fd3)
+
+    # check that for eps=0 power_eps equals standard np.power for N = 3.1
+
+    numba_f4, numba_fd4 = rlplib.power_eps( sv1, 3.1, eps)
+    numpy_f4, numpy_fd4 = np.power( sv1, 3.1), np.power( sv1, 3.1-1.0 )*3.1
+
+    assert np.allclose(numba_f4, numpy_f4)
+    assert np.allclose(numba_fd4, numpy_fd4)
+
+def test_power_eps5():
+
+    eps = 1e-8
+
+
+    sv1 = np.logspace(-6, 0, 51)
+
+    # check that N is clipped to zero when N < 0
+
+    numba_f1, numba_fd1 = rlplib.power_eps( sv1, -3.0, eps)
+    ref_f1, ref_fd1 = rlplib.power_eps( sv1, 0.0, eps)
+
+    assert np.allclose( numba_f1, ref_f1 )
+    assert np.allclose( numba_fd1, ref_fd1 )
+
+    # check that s is clipped to zero when s < 0, N = 3.0
+    sv2 = -np.logspace(-6, 0, 51)
+
+    numba_f2, numba_fd2 = rlplib.power_eps( sv2, 3.0, eps)
+    ref_f2, ref_fd2 = rlplib.power_eps( 0.0, 3.0, eps)
+
+    assert np.allclose( numba_f2, ref_f2 )
+    assert np.allclose( numba_fd2, ref_fd2 )
+
+    # check that s is clipped to zero when s < 0, N = 1.0
+    sv3 = -np.logspace(-6, 0, 51)
+
+    numba_f3, numba_fd3 = rlplib.power_eps( sv3, 1.0, eps)
+    ref_f3, ref_fd3 = rlplib.power_eps( 0.0, 1.0, eps)
+
+    assert np.allclose( numba_f3, ref_f3 )
+    assert np.allclose( numba_fd3, ref_fd3 )
+
+    # check that s is clipped to zero when s < 0, N = 0.5
+    sv4 = -np.logspace(-6, 0, 51)
+
+    numba_f4, numba_fd4 = rlplib.power_eps( sv4, 0.5, eps)
+    ref_f4, ref_fd4 = rlplib.power_eps( 0.0, 0.5, eps)
+
+    assert np.allclose( numba_f4, ref_f4 )
+    assert np.allclose( numba_fd4, ref_fd4 )
+
+    # check that s is clipped to zero when s < 0, N = 0.0
+    sv5 = -np.logspace(-6, 0, 51)
+
+    numba_f5, numba_fd5 = rlplib.power_eps( sv5, 0.0, eps)
+    ref_f5, ref_fd5 = rlplib.power_eps( 0.0, 0.0, eps)
+
+    assert np.allclose( numba_f5, ref_f5 )
+    assert np.allclose( numba_fd5, ref_fd5 )
+
+    # check that s is clipped to zero when s < 0, and N to zero when N = -1.0
+    sv6 = -np.logspace(-6, 0, 51)
+
+    numba_f6, numba_fd6 = rlplib.power_eps( sv6, -1.0, eps)
+    ref_f6, ref_fd6 = rlplib.power_eps( 0.0, 0.0, eps)
+
+    assert np.allclose( numba_f6, ref_f6 )
+    assert np.allclose( numba_fd6, ref_fd6 )

python/scallib001/tests/test_relperm_Corey1.py

@@ -0,0 +1,138 @@
+#------------------------------------------------------------------------------------------------
+#  Copyright (c) Shell Global Solutions International B.V. All rights reserved.
+#  Licensed under the MIT License. See License.txt in the project root for license information.
+#------------------------------------------------------------------------------------------------
+
+import pytest
+import numpy as np
+from scallib001 import relpermlib001 as rlplib
+
+
+def test_Corey_regular_derivatives1():
+
+    # When corey exponent < 1 infinite saturation derivatives occur at the end points, unless regulated using eps>0
+    swi = 0.20
+    sor = 0.15
+    krw = 0.35
+    kro = 0.80
+    lw = 0.5 # lw < 1, would cause infinite derivative at endpoint
+    lo = 0.5 # lo < 1, would cause infinite derivative at endpoint
+
+    eps = 1e-4
+
+    rlp_model = rlplib.Rlp2PCorey(
+        swi,
+        sor,
+        lw,
+        lo,
+        krw,
+        kro,
+        eps=eps,
+    )
+
+    # Choose to evaluate at the endpoints
+    swv = np.array([swi, 1-sor])
+
+    rlpw, rlpo, drlpw, drlpo = rlp_model.calc( swv )  
+
+    assert np.isfinite(drlpo[0])
+    assert np.isfinite(drlpw[1])
+
+def test_Corey_linear_model():
+
+    # Check derivatives for linear relperm
+    swi = 0.00
+    sor = 0.00
+    krw = 1.00
+    kro = 1.00
+    lw = 1.0
+    lo = 1.0
+
+    eps = 1e-4 
+
+    rlp_model = rlplib.Rlp2PCorey(
+        swi,
+        sor,
+        lw,
+        lo,
+        krw,
+        kro,
+        eps=eps,
+    )
+
+    # Choose to evaluate including the endpoints
+    swv = np.linspace(swi, 1-sor, 11)
+
+    rlpw, rlpo, drlpw, drlpo = rlp_model.calc( swv )  
+
+    print('rlpw', rlpw)
+    print('rlpo', rlpo)
+    print('drlpw', drlpw)
+    print('drlpo', drlpo)
+
+    assert np.allclose( rlpw, swv )
+    assert np.allclose( rlpo, 1.0-swv )
+    assert np.allclose( drlpw, +1.0 )
+    assert np.allclose( drlpo, -1.0 )
+
+def test_Corey_endpoints():
+
+    swi = 0.20
+    sor = 0.15
+    krw = 0.35
+    kro = 0.80
+    lw = 2.0
+    lo = 3.0
+
+    eps = 1e-4 
+
+    rlp_model = rlplib.Rlp2PCorey(
+        swi,
+        sor,
+        lw,
+        lo,
+        krw,
+        kro,
+        eps=eps,
+    )
+
+    # Choose to evaluate including the endpoints
+    swv = np.array([0.0, swi, 1-sor, 1.0])
+
+    rlpw, rlpo, drlpw, drlpo = rlp_model.calc( swv )  
+
+    assert np.allclose( rlpw, [0.0, 0.0, krw, krw] )
+    assert np.allclose( rlpo, [kro, kro, 0.0, 0.0] )
+    assert np.allclose( [drlpw[0], drlpw[-1]],  [0.0, 0.0] )
+    assert np.allclose( [drlpo[0], drlpo[-1]],  [0.0, 0.0] )
+
+def test_Corey_derivative1():
+
+    Sr1 = 0.20
+    Sr2 = 0.15
+    Ke1 = 0.35
+    Ke2 = 0.80
+    N1 = 2.4
+    N2 = 3.0
+
+    eps = 1e-4
+
+    rlp_model = rlplib.Rlp2PCorey(Sr1=Sr1,Sr2=Sr2,
+                                N1=N1,N2=N2,
+                                Ke1=Ke1,Ke2=Ke2,
+                                eps=eps)
+
+    deps = 0.01
+
+    swv = np.linspace(Sr1+deps, 1-Sr2-deps, 1001)
+
+    krw, kro, dkrw, dkro = rlp_model.calc(swv)
+
+    dkrw_num = np.gradient(krw, swv, edge_order=2)
+    dkro_num = np.gradient(kro, swv, edge_order=2)
+
+    print('dkrw abs error', np.abs(dkrw - dkrw_num).max(), 'rel error', np.abs((dkrw - dkrw_num)/dkrw).max())
+    print('dkro abs error', np.abs(dkro - dkro_num).max(), 'rel error', np.abs((dkro - dkro_num)/dkro).max())
+
+    assert np.allclose(dkrw, dkrw_num, atol=1e-3, rtol=1e-2)
+    assert np.allclose(dkro, dkro_num, atol=1e-3, rtol=1e-2)