From a7a7e3a23788424284c9693e67f78e001174d02d Mon Sep 17 00:00:00 2001 From: Vegard Gjeldvik Jervell Date: Wed, 17 Apr 2024 19:07:47 +0200 Subject: [PATCH 1/2] Doing division before applying exponent in potential functions. --- cpp/CMakeLists.txt | 2 +- cpp/MieKinGas.h | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index 29fb441..f1cfbb3 100644 --- a/cpp/CMakeLists.txt +++ b/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}") diff --git a/cpp/MieKinGas.h b/cpp/MieKinGas.h index 95dcae6..ba40542 100644 --- a/cpp/MieKinGas.h +++ b/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))); } From 3daadcf61cd5a77983f9a69bcd055a151e72fc28 Mon Sep 17 00:00:00 2001 From: Vegard Gjeldvik Jervell Date: Wed, 17 Apr 2024 19:34:16 +0200 Subject: [PATCH 2/2] Add test to ensure issue #25 remains solved. --- tests/README.md | 3 ++- tests/test_issues.py | 42 ++++++++++++++++++++++++++++++++++++++++++ tests/test_rdf.py | 2 +- 3 files changed, 45 insertions(+), 2 deletions(-) create mode 100644 tests/test_issues.py diff --git a/tests/README.md b/tests/README.md index 843c624..b901231 100644 --- a/tests/README.md +++ b/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. \ No newline at end of file +* `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. \ No newline at end of file diff --git a/tests/test_issues.py b/tests/test_issues.py new file mode 100644 index 0000000..809c529 --- /dev/null +++ b/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]) \ No newline at end of file diff --git a/tests/test_rdf.py b/tests/test_rdf.py index 318aca1..6ab6099 100644 --- a/tests/test_rdf.py +++ b/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) \ No newline at end of file + assert check_eq(rdf_tp, rdf_kgt, 1e-3) \ No newline at end of file