Merge pull request #26 from thermotools/much_repulse

Doing division before applying exponent in potential functions.

cpp/CMakeLists.txt

@@ -2,7 +2,7 @@ cmake_policy(VERSION 3.12)
 project(KineticGas LANGUAGES C CXX)
 
 cmake_minimum_required(VERSION 3.12)
-set(PYBIND11_PYTHON_VERSION "3.11")
+set(PYBIND11_PYTHON_VERSION "3")
 string(ASCII 27 Esc)
 message("${Esc}[34mPYBIND11_PYTHON_VERSION is : ${PYBIND11_PYTHON_VERSION}")
 message("${Esc}[34mPYTHON_EXECUTABLE is : ${PYTHON_EXECUTABLE}")

cpp/MieKinGas.h

@@ -46,13 +46,13 @@ class MieKinGas : public Spherical {
         return C[i][j] * eps[i][j] * (pow(sigma[i][j] / r, lr[i][j]) - pow(sigma[i][j] / r, la[i][j]));
     }
     double potential_derivative_r(int i, int j, double r) override {
-        return C[i][j] * eps[i][j] * ((la[i][j] * pow(sigma[i][j], la[i][j]) / pow(r, la[i][j] + 1)) 
-                                        - (lr[i][j] * pow(sigma[i][j], lr[i][j]) / pow(r, lr[i][j] + 1)));
+        return C[i][j] * eps[i][j] * ((la[i][j] * pow(sigma[i][j] / r, la[i][j]) / r)
+                                        - (lr[i][j] * pow(sigma[i][j] / r, lr[i][j]) / r));
     }
 
     double potential_dblderivative_rr(int i, int j, double r) override {
-        return C[i][j] * eps[i][j] * ((lr[i][j] * (lr[i][j] + 1) * pow(sigma[i][j], lr[i][j]) / pow(r, lr[i][j] + 2)) 
-                                    - (la[i][j] * (la[i][j] + 1) * pow(sigma[i][j], la[i][j]) / pow(r, la[i][j] + 2)));
+        return C[i][j] * eps[i][j] * ((lr[i][j] * (lr[i][j] + 1) * pow(sigma[i][j] / r, lr[i][j]) / pow(r, 2))
+                                    - (la[i][j] * (la[i][j] + 1) * pow(sigma[i][j] / r, la[i][j]) / pow(r, 2)));
     }
 
 

tests/README.md

@@ -42,4 +42,5 @@ Current test modules are:
 * `test_binary_component_order.py` - Checks that order of components in binary mixtures does not affect output
 * `test_summation_consistency.py` - Checks that various things that should sum to 1 or 0 do that.
 * `test_binary_limits.py` - Checks that ternary coefficients reduce to the appropriate binary coefficients in the binary limit.
-* `test_flux_transforms.py` - Checks that the frame of reference transformations give fluxes that appropriately sum to zero.
+* `test_flux_transforms.py` - Checks that the frame of reference transformations give fluxes that appropriately sum to zero.
+* `test_issues` - New test added anytime an issue is resolved, to ensure the issue remains resolved.

tests/test_issues.py

@@ -0,0 +1,42 @@
+from pykingas.MieKinGas import MieKinGas
+from scipy.constants import Avogadro, Boltzmann
+import numpy as np
+from tools import check_eq
+import pytest
+
+@pytest.mark.parametrize('lr', [20, 30, 40, 50])
+def test_very_repulse(lr):
+    """
+    Issue #25 (https://github.com/thermotools/KineticGas/issues/25)
+    Check that viscosity and thermal conductivity run for highly repulsive systems.
+    """
+    m = 10
+    mie = MieKinGas('LJF', mole_weights=[m, m], lr=[lr, lr])
+
+    sigma = mie.sigma[0][0]
+    eps = mie.epsilon[0]
+    eps_div_k = eps / Boltzmann
+
+    test_visc_vals = {20 : 0.2136582471368343,
+                        30 : 0.21469802292722762,
+                        40 : 0.21527111014341763,
+                        50 : 0.21563637353012097}
+
+    test_cond_vals = {20 : 0.7991821474987899,
+                        30 : 0.8022966936925214,
+                        40 : 0.8039439619090943,
+                        50 : 0.8049669080299001}
+
+    T_red = 2.0
+    rho_red = 0.1
+    T = T_red * eps_div_k
+    rho = rho_red * Avogadro * sigma**3
+
+    visc_unit = np.sqrt(eps * (m * 1e-3 / Avogadro)) / sigma**2
+    cond_unit = Boltzmann * np.sqrt(eps/ (m * 1e-3 / Avogadro)) / sigma**2
+
+    visc = mie.viscosity(T, 1 / rho, [0.5, 0.5], N=2) / visc_unit
+    cond = mie.thermal_conductivity(T, 1 / rho, [0.5, 0.5], N=2) / cond_unit
+
+    assert check_eq(visc, test_visc_vals[lr])
+    assert check_eq(cond, test_cond_vals[lr])

tests/test_rdf.py

@@ -37,4 +37,4 @@ def test_vs_thermopack(T, rho, lr):
         if (vl < 1 / rho) and (1 / rho < vg):
             return # Two-phase region
 
-    assert check_eq(rdf_tp, rdf_kgt, 1e-8)
+    assert check_eq(rdf_tp, rdf_kgt, 1e-3)