diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index 3a628c7..29fb441 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -1,7 +1,8 @@ +cmake_policy(VERSION 3.12) project(KineticGas LANGUAGES C CXX) cmake_minimum_required(VERSION 3.12) -# set(PYBIND11_PYTHON_VERSION "3.10") +set(PYBIND11_PYTHON_VERSION "3.11") string(ASCII 27 Esc) message("${Esc}[34mPYBIND11_PYTHON_VERSION is : ${PYBIND11_PYTHON_VERSION}") message("${Esc}[34mPYTHON_EXECUTABLE is : ${PYTHON_EXECUTABLE}") @@ -32,4 +33,4 @@ else() endif() add_subdirectory(${PYBIND11_ROOT} release) -pybind11_add_module(${TARGET} ${SOURCES}) +pybind11_add_module(${TARGET} ${SOURCES}) \ No newline at end of file diff --git a/cpp/HardSphere.h b/cpp/HardSphere.h index e79c97b..04ca469 100644 --- a/cpp/HardSphere.h +++ b/cpp/HardSphere.h @@ -18,13 +18,13 @@ class HardSphere : public KineticGas { : KineticGas(mole_weights, is_idealgas), sigma{sigmaij}{ } - double omega(const int& i, const int& j, const int& l, const int& r, const double& T) override { + double omega(int i, int j, int l, int r, double T) override { double w = w_integral(i, j, T, l, r); if (i == j) return pow(sigma.at(i).at(j), 2) * sqrt((PI * BOLTZMANN * T) / m.at(i)) * w; return 0.5 * pow(sigma.at(i).at(j), 2) * sqrt(2 * PI * BOLTZMANN * T / (m0[i][j] * M[i][j] * M[j][i])) * w; } - double w_integral(const int& i, const int& j, const double& T, const int& l, const int& r){ + double w_integral(int i, int j, double T, int l, int r){ int f = Fac(r + 1).eval(); if (l % 2 == 0){ return 0.25 * (2 - ((1.0 / (l + 1)) * 2)) * f; @@ -59,6 +59,6 @@ class HardSphere : public KineticGas { return rdf; } - std::vector> get_contact_diameters(double rho, double T, const std::vector& x) override {return sigma;} + std::vector> get_collision_diameters(double rho, double T, const std::vector& x) override {return sigma;} }; \ No newline at end of file diff --git a/cpp/KineticGas.cpp b/cpp/KineticGas.cpp index 3c4fc07..4202a21 100644 --- a/cpp/KineticGas.cpp +++ b/cpp/KineticGas.cpp @@ -63,7 +63,7 @@ std::vector KineticGas::get_K_factors(double rho, double T, const std::v std::vector> rdf = get_rdf(rho, T, mole_fracs); std::vector K(Ncomps, 0.); - std::vector> cd = get_contact_diameters(rho, T, mole_fracs); + std::vector> cd = get_collision_diameters(rho, T, mole_fracs); for (int i = 0; i < Ncomps; i++){ for (int j = 0; j < Ncomps; j++){ K[i] += mole_fracs[j] * pow(cd[i][j], 3) * M[i][j] * M[j][i] * rdf[i][j]; @@ -80,7 +80,7 @@ std::vector KineticGas::get_K_prime_factors(double rho, double T, const std::vector> rdf = get_rdf(rho, T, mole_fracs); std::vector K_prime(Ncomps, 0.); - std::vector> cd = get_contact_diameters(rho, T, mole_fracs); + std::vector> cd = get_collision_diameters(rho, T, mole_fracs); for (int i = 0; i < Ncomps; i++){ for (int j = 0; j < Ncomps; j++){ K_prime[i] += mole_fracs[j] * M[j][i] * pow(cd[i][j], 3.) * rdf[i][j]; @@ -293,7 +293,7 @@ std::vector KineticGas::get_viscosity_vector(double rho, double T, const // A_pqrl factors, used for conductivity and diffusion // // --------------------------------------------------------------------------------------------------- // -double KineticGas::A(const int& p, const int& q, const int& r, const int& l){ +double KineticGas::A(int p, int q, int r, int l) const { double value{0.0}; int max_i = std::min(std::min(p, q), std::min(r, p + q + 1 - r)); for (int i = l - 1; i <= max_i; i++){ @@ -304,8 +304,7 @@ double KineticGas::A(const int& p, const int& q, const int& r, const int& l){ return value; } -double KineticGas::A_prime(const int& p, const int& q, const int& r, const int& l, - const double& tmp_M1, const double& tmp_M2){ +double KineticGas::A_prime(int p, int q, int r, int l, double tmp_M1, double tmp_M2) const { double F = (pow(tmp_M1, 2) + pow(tmp_M2, 2)) / (2 * tmp_M1 * tmp_M2); double G = (tmp_M1 - tmp_M2) / tmp_M2; @@ -338,12 +337,11 @@ double KineticGas::A_prime(const int& p, const int& q, const int& r, const int& // B_pqrl factors, used for viscosity // // --------------------------------------------------------------------------------------------------- // -inline Frac poch(const int& z, const int& n){ // Pochhammer notation (z)_n, used in the B-factors +inline Frac poch(int z, int n){ // Pochhammer notation (z)_n, used in the B-factors return Fac(z + n - 1) / Fac(z - 1); } -double KineticGas::B_prime(const int& p, const int& q, const int& r, const int& l, - const double& M1, const double& M2){ +double KineticGas::B_prime(int p, int q, int r, int l, double M1, double M2) const { if (r == p + q + 3) return 0.0; double val{0.0}; double inner{0.0}; @@ -381,8 +379,7 @@ double KineticGas::B_prime(const int& p, const int& q, const int& r, const int& return val; } -double KineticGas::B_dblprime(const int& p, const int& q, const int& r, const int& l, - const double& M1, const double& M2){ +double KineticGas::B_dblprime(int p, int q, int r, int l, double M1, double M2) const { if (r == p + q + 3) return 0.0; double val{0.0}; Frac num; @@ -420,7 +417,7 @@ double KineticGas::B_dblprime(const int& p, const int& q, const int& r, const in // H-integrals, used for conductivity and diffusion // // --------------------------------------------------------------------------------------------------- // -double KineticGas::H_ij(const int& p, const int& q, const int& i, const int& j, const double& T){ +double KineticGas::H_ij(int p, int q, int i, int j, double T){ double tmp_M1{M[i][j]}, tmp_M2{M[j][i]}; double value{0.0}; int max_l = std::min(p, q) + 1; @@ -435,7 +432,7 @@ double KineticGas::H_ij(const int& p, const int& q, const int& i, const int& j, return value; } -double KineticGas::H_i(const int& p, const int& q, const int& i, const int& j, const double& T){ +double KineticGas::H_i(int p, int q, int i, int j, double T){ double tmp_M1{M[i][j]}, tmp_M2{M[j][i]}; double value{0.0}; @@ -455,7 +452,7 @@ double KineticGas::H_i(const int& p, const int& q, const int& i, const int& j, c // L-integrals, used for viscosity // // --------------------------------------------------------------------------------------------------- // -double KineticGas::L_ij(const int& p, const int& q, const int& i, const int& j, const double& T){ +double KineticGas::L_ij(int p, int q, int i, int j, double T){ double val{0.0}; double M1{M[i][j]}; double M2{M[j][i]}; @@ -467,7 +464,7 @@ double KineticGas::L_ij(const int& p, const int& q, const int& i, const int& j, val *= 16.0 / 3.0; return val; } -double KineticGas::L_i(const int& p, const int& q, const int& i, const int& j, const double& T){ +double KineticGas::L_i(int p, int q, int i, int j, double T){ double val{0.0}, M1{M[i][j]}, M2{M[j][i]}; for (int l = 1; l <= std::min(p, q) + 2; l++){ for (int r = l; r <= p + q + 4 - l; r++){ diff --git a/cpp/KineticGas.h b/cpp/KineticGas.h index 01cab10..1d51585 100644 --- a/cpp/KineticGas.h +++ b/cpp/KineticGas.h @@ -64,62 +64,34 @@ struct OmegaPoint{ class KineticGas{ public: - const unsigned long Ncomps; // must be ulong to be initialised from mole_weights.size() - const bool is_idealgas; - std::vector m; - std::vector> M, m0; - std::map omega_map; KineticGas(std::vector mole_weights, bool is_idealgas); virtual ~KineticGas(){}; // Collision integrals - virtual double omega(const int& i, const int& j, const int& l, const int& r, const double& T) = 0; + virtual double omega(int i, int j, int l, int r, double T) = 0; // The "distance between particles" at contact. // NOTE: Will Return [0, 0, ... 0] for models with is_idealgas=True. - virtual std::vector> get_contact_diameters(double rho, double T, const std::vector& x) = 0; - - /* - The radial distribution function "at contact" for the given potential model - model_rdf is only called if object is initialized with is_idealgas=true - If a potential model is implemented only for the ideal gas state, its implementation - of model_rdf should throw an std::invalid_argument error. - */ - virtual std::vector> model_rdf(double T, double rho, const std::vector& mole_fracs) = 0; + virtual std::vector> get_collision_diameters(double rho, double T, const std::vector& x) = 0; - std::vector get_wt_fracs(const std::vector mole_fracs); // Compute weight fractions from mole fractions + // Radial distribution function "at contact". Inheriting classes must implement model_rdf. std::vector> get_rdf(double rho, double T, const std::vector& mole_fracs) { if (is_idealgas) return std::vector>(Ncomps, std::vector(Ncomps, 1.)); return model_rdf(rho, T, mole_fracs); } + /* + The radial distribution function "at contact" for the given potential model. model_rdf is only called if object + is initialized with is_idealgas=true. If a potential model is implemented only for the ideal gas state, its + implementation of model_rdf should throw an std::invalid_argument error. + */ + virtual std::vector> model_rdf(double rho, double T, const std::vector& mole_fracs) = 0; -// ------------------------------------------------------------------------------------------------------------------------ // -// -------------------------------------------------- Square bracket integrals -------------------------------------------- // - - // Linear combination weights by Tompson, Tipton and Lloyalka - double A(const int& p, const int& q, const int& r, const int& l); - double A_prime(const int& p, const int& q, const int& r, const int& l, const double& tmp_M1, const double& tmp_M2); - double A_trippleprime(const int& p, const int& q, const int& r, const int& l); - - // The diffusion and conductivity related square bracket integrals - double H_ij(const int& p, const int& q, const int& i, const int& j, const double& T); // [S^(p)_{3/2}(U^2_i), S^(q)_{3/2}(U^2_j)]_{ij} - double H_i(const int& p, const int& q, const int& i, const int& j, const double& T); // [S^(p)_{3/2}(U^2_i), S^(q)_{3/2}(U^2_i)]_{ij} - double H_simple(const int& p, const int& q, const int& i, const double& T); // [S^(p)_{3/2}(U^2_i), S^(q)_{3/2}(U^2_i)]_{i} - - // Linear combination weights by Tompson, Tipton and Lloyalka - double B_prime(const int& p, const int& q, const int& r, const int& l, - const double& M1, const double& M2); - double B_dblprime(const int& p, const int& q, const int& r, const int& l, - const double& M1, const double& M2); - double B_trippleprime(const int& p, const int& q, const int& r, const int& l); - - // Viscosity related square bracket integrals - double L_ij(const int& p, const int& q, const int& i, const int& j, const double& T); // [S^(p)_{5/2}(U^2_i), S^(q)_{5/2}(U^2_j)]_{ij} - double L_i(const int& p, const int& q, const int& i, const int& j, const double& T); // [S^(p)_{5/2}(U^2_i), S^(q)_{5/2}(U^2_i)]_{ij} - double L_simple(const int& p, const int& q, const int& i, const double& T); // [S^(p)_{5/2}(U^2_i), S^(q)_{5/2}(U^2_i)]_{i} + std::vector get_wt_fracs(const std::vector mole_fracs); // Compute weight fractions from mole fractions -// ---------------------------------------------------------------------------------------------------------------------------------------------- // -// ---------------------------------- Matrices and vectors for multicomponent, density corrected solutions -------------------------------------- // + // ------------------------------------------------------------------------------------------------------------------------- // + // ----- Matrices and vectors for the sets of equations (6-10) in Revised Enskog Theory for Mie fluids -------------------- // + // ----- doi : 10.1063/5.0149865, which are solved to obtain the Sonine polynomial expansion coefficients ------------------ // + // ----- for the velocity distribution functions. -------------------------------------------------------------------------- // std::vector> get_conductivity_matrix(double rho, double T, const std::vector& x, int N); std::vector get_diffusion_vector(double rho, double T, const std::vector& x, int N); @@ -128,32 +100,70 @@ class KineticGas{ std::vector get_conductivity_vector(double rho, double T, const std::vector& x, int N); std::vector> get_viscosity_matrix(double rho, double T, const std::vector&x, int N); std::vector get_viscosity_vector(double rho, double T, const std::vector& x, int N); - - // Eq. (1.2) of 'multicomponent docs' - std::vector get_K_factors(double rho, double T, const std::vector& mole_fracs); - // Eq. (5.4) of 'multicomponent docs' - std::vector get_K_prime_factors(double rho, double T, const std::vector& mole_fracs); + std::vector get_K_factors(double rho, double T, const std::vector& mole_fracs); // Eq. (1.2) of 'multicomponent docs' + std::vector get_K_prime_factors(double rho, double T, const std::vector& mole_fracs); // Eq. (5.4) of 'multicomponent docs' + +// ------------------------------------------------------------------------------------------------------------------------ // +// --------------------------------------- KineticGas internals are below here -------------------------------------------- // +// -------------------------------- End users should not need to care about any of this ----------------------------------- // + + protected: + const size_t Ncomps; + const bool is_idealgas; + std::vector m; + std::vector> M, m0; + std::map omega_map; // ----------------------------------------------------------------------------------------------------------------------------------- // // --------------------------------------- Methods to facilitate multithreading ------------------------------------------------------ // /* - Filling the A-matrix is the most computationally intensive part of the module - The matrix elements may be computed independently. Therefore the functions fill_A_matrix_() have been - implemented to facilitate multithreading. j indicates the number of functions the matrix generation has been split - over, i differentiates between functions. - So fill_A_matrix_i4 is a set of four functions intended to be run on four separate threads, while - fill_A_matrix_i2 is a set of two functions intended to be run on two separate threads. - For a graphical representation of how the functions split the matrix, see the implementation file. + Computing collision integrals is the most computationally intensive part of computing a transport property, + given that one does not use correlations for the computation. Therefore, computed collision integrals are stored + in this->omega_map. These methods are used to deduce which collision integrals are required to compute a given + property, and then split the computation of those integrals among several threads. When computation of the + property proceeds, the integrals are then retrieved from this->omega_map. + + To adjust the number of threads used to compute collision integrals, set the compile-time constant Ncores in + KineticGas_mthr.cpp. In practice, very little is to be gained by increasing this beyond 10, unless you are + using high Enskog approximation orders (>4), or are working with a large number of components (because there + is no point in using more threads than required integrals). + + These need to be protected instead of private, because QuantumMie needs to override them to prevent a race condition + without locking everything with a mutex on every call to omega. */ - void precompute_omega(const std::vector& i_vec, const std::vector& j_vec, + virtual void precompute_conductivity_omega(int N, double T); + virtual void precompute_viscosity_omega(int N, double T); + virtual void precompute_omega(const std::vector& i_vec, const std::vector& j_vec, const std::vector& l_vec, const std::vector& r_vec, double T); - void precompute_conductivity_omega(int N, double T); - void precompute_diffusion_omega(int N, double T); - void precompute_th_diffusion_omega(int N, double T); - void precompute_viscosity_omega(int N, double T); + private: + void precompute_diffusion_omega(int N, double T); // Forwards call to precompute_conductivity_omega. Override that instead. + void precompute_th_diffusion_omega(int N, double T); // Forwards call to precompute_conductivity_omega. Override that instead. + private: +// ------------------------------------------------------------------------------------------------------------------------ // +// ---------------------------------------------- Square bracket integrals ------------------------------------------------ // + + // Linear combination weights by Tompson, Tipton and Lloyalka + double A(int p, int q, int r, int l) const; + double A_prime(int p, int q, int r, int l, double tmp_M1, double tmp_M2) const; + double A_trippleprime(int p, int q, int r, int l) const; + + // The diffusion and conductivity related square bracket integrals + double H_ij(int p, int q, int i, int j, double T); // [S^(p)_{3/2}(U^2_i), S^(q)_{3/2}(U^2_j)]_{ij} + double H_i(int p, int q, int i, int j, double T); // [S^(p)_{3/2}(U^2_i), S^(q)_{3/2}(U^2_i)]_{ij} + double H_simple(int p, int q, int i, double T); // [S^(p)_{3/2}(U^2_i), S^(q)_{3/2}(U^2_i)]_{i} + + // Linear combination weights by Tompson, Tipton and Lloyalka + double B_prime(int p, int q, int r, int l, double M1, double M2) const; + double B_dblprime(int p, int q, int r, int l, double M1, double M2) const; + double B_trippleprime(int p, int q, int r, int l) const; + + // Viscosity related square bracket integrals + double L_ij(int p, int q, int i, int j, double T); // [S^(p)_{5/2}(U^2_i), S^(q)_{5/2}(U^2_j)]_{ij} + double L_i(int p, int q, int i, int j, double T); // [S^(p)_{5/2}(U^2_i), S^(q)_{5/2}(U^2_i)]_{ij} + double L_simple(int p, int q, int i, double T); // [S^(p)_{5/2}(U^2_i), S^(q)_{5/2}(U^2_i)]_{i} }; inline int delta(int i, int j) {return (i == j) ? 1 : 0;} \ No newline at end of file diff --git a/cpp/MieKinGas.cpp b/cpp/MieKinGas.cpp index 117d6f9..24fa5cd 100644 --- a/cpp/MieKinGas.cpp +++ b/cpp/MieKinGas.cpp @@ -8,8 +8,8 @@ See : J. Chem. Phys. 139, 154504 (2013); https://doi.org/10.1063/1.4819786 using namespace mie_rdf_constants; -double MieKinGas::omega(const int& i, const int& j, const int& l, const int& r, const double& T){ - if ((l <= 2) && (r <= 3) && (abs(la[i][j] - 6.) < 1e-10)){ // Use Correlation by Fokin, Popov and Kalashnikov, High Temperature, Vol. 37, No. 1 (1999) +double MieKinGas::omega(int i, int j, int l, int r, double T){ + if ((l <= 2) && (r >= l) && (r <= 3) && (abs(la[i][j] - 6.) < 1e-10)){ // Use Correlation by Fokin, Popov and Kalashnikov, High Temperature, Vol. 37, No. 1 (1999) // NOTE: There is a typo in Eq. (4b) in the article, ln(1/m) should be ln(m). // See: I. H. Bell et. al. J. Chem. Eng. Data (2020), https://doi.org/10.1021/acs.jced.9b00455 // The correlation implemented here has been verified vs. tabulated data. @@ -64,14 +64,14 @@ double MieKinGas::omega_recursive_factor(int i, int j, int l, int r, double T_st return 1. + dlnomega_dlnT / (r + 2); } -std::vector> MieKinGas::model_rdf(double rho, double T, const std::vector& x){ +std::vector> MieKinGas::saft_rdf(double rho, double T, const std::vector& x, int order, bool g2_correction){ double beta = (1. / (BOLTZMANN * T)); std::vector> d_BH = get_BH_diameters(T); std::vector> x0 = get_x0(d_BH); - std::vector> gHS = rdf_HS(rho, T, x); - std::vector> g1 = rdf_g1_func(rho, x, d_BH); - std::vector> g2 = rdf_g2_func(rho, T, x, d_BH, x0); - std::vector> g(Ncomps, std::vector(Ncomps)); + std::vector> gHS = rdf_HS(rho, x, d_BH); + std::vector> g1 = (order > 0) ? rdf_g1_func(rho, x, d_BH) : std::vector>(Ncomps, std::vector(Ncomps, 0.)); + std::vector> g2 = (order > 1) ? rdf_g2_func(rho, T, x, d_BH, x0, g2_correction) : std::vector>(Ncomps, std::vector(Ncomps, 0.)); + std::vector> g(Ncomps, std::vector(Ncomps, 0.0)); for (int i = 0; i < Ncomps; i++){ for (int j = i; j < Ncomps; j++){ g[i][j] = gHS[i][j] * exp(beta * eps[i][j] * (g1[i][j] / gHS[i][j]) + pow(beta * eps[i][j], 2) * (g2[i][j] / gHS[i][j])); @@ -80,9 +80,7 @@ std::vector> MieKinGas::model_rdf(double rho, double T, cons } return g; } - -std::vector> MieKinGas::rdf_HS(double rho, double T, const std::vector& x){ - std::vector> d_BH = get_BH_diameters(T); +std::vector> MieKinGas::rdf_HS(double rho, const std::vector& x, const std::vector>& d_BH){ double zeta = zeta_x_func(rho, x, d_BH); double k0 = - log(1 - zeta) + (42. * zeta - 39. * pow(zeta, 2) + 9. * pow(zeta, 3) - 2 * pow(zeta, 4)) / (6 * pow(1 - zeta, 3)); double k1 = (pow(zeta, 4) + 6. * pow(zeta, 2) - 12. * zeta) / (2. * pow(1 - zeta, 3)); @@ -102,6 +100,11 @@ std::vector> MieKinGas::rdf_HS(double rho, double T, const s return rdf; } +std::vector> MieKinGas::rdf_HS(double rho, double T, const std::vector& x){ + std::vector> d_BH = get_BH_diameters(T); + return rdf_HS(rho, x, d_BH); +} + std::vector> MieKinGas::rdf_g1_func(double rho, const std::vector& x, const std::vector>& d_BH) { @@ -130,7 +133,8 @@ std::vector> MieKinGas::rdf_g1_func(double rho, double T, co std::vector> MieKinGas::rdf_g2_func(double rho, double T, const std::vector& x, const std::vector>& d_BH, - const std::vector>& x0) + const std::vector>& x0, + bool g2_correction) { std::vector> g2(Ncomps, std::vector(Ncomps)); std::vector> _la(Ncomps, std::vector(Ncomps)); @@ -158,7 +162,7 @@ std::vector> MieKinGas::rdf_g2_func(double rho, double T, co std::vector> da2ij_div_chi_drho = da2ij_div_chi_drho_func(rho, x, K_HS, d_BH, x0); std::vector> a2ij = a2ij_func(rho, x, K_HS, rdf_chi_HS, d_BH, x0); double zeta_x_HS = zeta_x_func(rho, x, sigma); - std::vector> gamma_c = gamma_corr(zeta_x_HS, T); + std::vector> gamma_c = g2_correction ? gamma_corr(zeta_x_HS, T) : std::vector>(Ncomps, std::vector(Ncomps, 0.)); for (int i = 0; i < Ncomps; i++){ for (int j = i; j < Ncomps; j++){ g2[i][j] = (1 + gamma_c[i][j]) @@ -176,10 +180,10 @@ std::vector> MieKinGas::rdf_g2_func(double rho, double T, co return g2; } -std::vector> MieKinGas::rdf_g2_func(double rho, double T, const std::vector& x){ +std::vector> MieKinGas::rdf_g2_func(double rho, double T, const std::vector& x, bool g2_correction){ std::vector> d_BH = get_BH_diameters(T); std::vector> x0 = get_x0(d_BH); - return rdf_g2_func(rho, T, x, d_BH, x0); + return rdf_g2_func(rho, T, x, d_BH, x0, g2_correction); } std::vector> MieKinGas::get_b_max(double T){ @@ -219,7 +223,7 @@ std::vector> MieKinGas::get_b_max(double T){ return bmax; } -std::vector> MieKinGas::get_contact_diameters(double rho, double T, const std::vector& x){ +std::vector> MieKinGas::get_collision_diameters(double rho, double T, const std::vector& x){ // Evaluate the integral of Eq. (40) in RET for Mie fluids (doi: 10.1063/5.0149865) // Note: For models with is_idealgas==true, zeros are returned, as there is no excluded volume at infinite dilution. // Weights and nodes for 6-point Gauss-Legendre quadrature @@ -249,14 +253,15 @@ std::vector> MieKinGas::get_contact_diameters(double rho, do } std::vector> MieKinGas::get_BH_diameters(double T){ - // Gauss-Legendre points taken from SAFT-VR-MIE docs (see: ThermoPack) + // 20-point Gauss-Legendre from 0.5 sigma to 1 sigma. Using constant value of 1 for integrand within (0, 0.5) sigma. std::vector> d_BH(Ncomps, std::vector(Ncomps, 0.0)); double beta = 1. / (BOLTZMANN * T); for (int i = 0; i < Ncomps; i++){ - for (int n = 0; n < 10; n++){ - d_BH[i][i] += gl_w[n] * (1. - exp(- beta * potential(i, i, sigma[i][i] * (gl_x[n] + 1) / 2.))); + d_BH[i][i] = 2.; + for (int n = 0; n < 20; n++){ + d_BH[i][i] += gl_w[n] * (1. - exp(- beta * potential(i, i, sigma[i][i] * (gl_x[n] / 4. + 3. / 4.)))); } - d_BH[i][i] *= sigma[i][i] / 2; + d_BH[i][i] *= sigma[i][i] / 4.; } for (int i = 0; i < Ncomps - 1; i++){ for (int j = i + 1; j < Ncomps; j++){ @@ -288,7 +293,8 @@ std::vector> MieKinGas::a_1s_func(double rho, for (int i = 0; i < Ncomps; i++){ for (int j = i; j < Ncomps; j++){ zeta_eff = zeta_eff_func(rho, x, d_BH, lambda[i][j]); - a1s[i][j] = (-2. * rho * PI * eps[i][j] * pow(d_BH[i][j], 3) / (lambda[i][j] - 3)) * (1. - zeta_eff/2.) / pow(1 - zeta_eff, 3); + a1s[i][j] = (-2. * rho * PI * eps[i][j] * pow(d_BH[i][j], 3) / (lambda[i][j] - 3)) + * (1. - zeta_eff/2.) / pow(1. - zeta_eff, 3); a1s[j][i] = a1s[i][j]; } } @@ -375,6 +381,13 @@ std::vector> MieKinGas::B_func(double rho, const std::vector return B; } +std::vector> MieKinGas::B_func(double rho, double T, const std::vector& x, + const std::vector>& lambda) +{ + std::vector> d_BH = get_BH_diameters(T); + return B_func(rho, x, d_BH, lambda); +} + std::vector> MieKinGas::dBdrho_func(double rho, const std::vector& x, const std::vector>& d_BH, const std::vector>& lambda) @@ -388,14 +401,23 @@ std::vector> MieKinGas::dBdrho_func(double rho, const std::v for (int j = i; j < Ncomps; j++){ double dzxdrho = dzetax_drho_func(x, d_BH); double zx = zeta_x_func(rho, x, d_BH); - dBdrho[i][j] = B[i][j] / rho + 2. * PI * eps[i][j] * d_BH[i][j] * dzxdrho - * (I[i][j] * (2.5 - zx) / pow(1 - zx, 4) - - 4.5 * J[i][j] * (1 + 3 * zx - pow(zx, 2)) / pow(1 - zx, 2)); + dBdrho[i][j] = B[i][j] / rho + + 2. * PI * rho * eps[i][j] * pow(d_BH[i][j], 3) + * (I[i][j] * (- 0.5 * (1 - zx) + 3. * (1. - 0.5 * zx)) / pow(1 - zx, 4) + - J[i][j] * (9. * (1. + 2. * zx) * (1. - zx) + 27. * zx * (1. + zx)) / (2. * pow(1 - zx, 4)) + ) * dzxdrho; } } return dBdrho; } +std::vector> MieKinGas::dBdrho_func(double rho, double T, const std::vector& x, + const std::vector>& lambda) +{ + std::vector> d_BH = get_BH_diameters(T); + return dBdrho_func(rho, x, d_BH, lambda); +} + std::vector> MieKinGas::a1ij_func(double rho, double T, const std::vector& x) { std::vector> d_BH = get_BH_diameters(T); @@ -613,9 +635,11 @@ std::vector> MieKinGas::gamma_corr(double zeta_x, double T){ for (int i = 0; i < Ncomps; i++){ for (int j = i; j < Ncomps; j++){ - double theta = exp(eps[i][j] / (BOLTZMANN * T)) + 1.; + double theta = exp(eps[i][j] / (BOLTZMANN * T)) - 1.; gamma[i][j] = phi[0] * (1 - tanh(phi[1] * (phi[2] - alpha[i][j]))) * zeta_x * theta * exp(phi[3] * zeta_x + phi[4] * pow(zeta_x, 2)); + gamma[j][i] = gamma[i][j]; + } } return gamma; @@ -626,8 +650,8 @@ double MieKinGas::zeta_x_func(double rho, const std::vector>& d_BH) { double zeta{0.0}; - for (int i = 0; i < x.size(); i++){ - for (int j = 0; j < x.size(); j++){ + for (int i = 0; i < Ncomps; i++){ + for (int j = 0; j < Ncomps; j++){ zeta += x[i] * x[j] * pow(d_BH[i][j], 3); } } diff --git a/cpp/MieKinGas.h b/cpp/MieKinGas.h index b015f57..95dcae6 100644 --- a/cpp/MieKinGas.h +++ b/cpp/MieKinGas.h @@ -56,13 +56,13 @@ class MieKinGas : public Spherical { } - double omega(const int& i, const int& j, const int& l, const int& r, const double& T) override; + double omega(int i, int j, int l, int r, double T) override; double omega_correlation(int i, int j, int l, int r, double T_star); double omega_recursive_factor(int i, int j, int l, int r, double T); // The hard sphere integrals are used as the reducing factor for the correlations. // So we need to compute the hard-sphere integrals to convert the reduced collision integrals from the // Correlation by Fokin et. al. to the "real" collision integrals. - inline double omega_hs(const int& i, const int& j, const int& l, const int& r, const double& T){ + inline double omega_hs(int i, int j, int l, int r, double T){ double w = PI * pow(sigma[i][j], 2) * 0.5 * (r + 1); for (int ri = r; ri > 1; ri--) {w *= ri;} if (l % 2 == 0){ @@ -96,19 +96,35 @@ class MieKinGas : public Spherical { // bmax[i][j] is in units of sigma[i][j] // bmax = The maximum value of the impact parameter at which deflection angle (chi) is positive std::vector> get_b_max(double T); - std::vector> get_contact_diameters(double rho, double T, const std::vector& x) override; + std::vector> get_collision_diameters(double rho, double T, const std::vector& x) override; virtual std::vector> get_BH_diameters(double T); + std::vector> get_vdw_alpha(){return alpha;} // Methods for computing the radial distribution function at contact - std::vector> model_rdf(double rho, double T, const std::vector& x) override; + // Note: A lot of these methods have two overloads: One that takes the temperature and density, and comptes + // the BH diameter, packing fractions etc, and another that takes the BH diameters, packing fractions, + // and other variables that must be pre-computed directly. I'm not sure which version of these is in primary use + // In the "standard" call chain when `model_rdf` is called, but someone should at some point ensure that we + // are not computing a bunch of unnessecary BH diameters. + // Note: For SAFT-type models, we + inline std::vector> model_rdf(double rho, double T, const std::vector& x) override { + return saft_rdf(rho, T, x, 2); + } + + // To directly compute the RDF at different pertubation orders. Not used in property computations. + std::vector> saft_rdf(double rho, double T, const std::vector& x, int order=2, bool g2_correction=true); + + std::vector> rdf_HS(double rho, const std::vector& x, + const std::vector>& d_BH); std::vector> rdf_HS(double rho, double T, const std::vector& x); std::vector> rdf_g1_func(double rho, const std::vector& x, const std::vector>& d_BH); std::vector> rdf_g1_func(double rho, double T, const std::vector& x); std::vector> rdf_g2_func(double rho, double T, const std::vector& x, const std::vector>& d_BH, - const std::vector>& x0); - std::vector> rdf_g2_func(double rho, double T, const std::vector& x); + const std::vector>& x0, + bool g2_correction=true); + std::vector> rdf_g2_func(double rho, double T, const std::vector& x, bool g2_correction=true); std::vector> get_x0(const std::vector>& d_BH); std::vector> a_1s_func(double rho, const std::vector& x, const std::vector>& d_BH, @@ -128,9 +144,13 @@ class MieKinGas : public Spherical { std::vector> B_func(double rho, const std::vector& x, const std::vector>& d_BH, const std::vector>& lambda); + std::vector> B_func(double rho, double T, const std::vector& x, + const std::vector>& lambda); std::vector> dBdrho_func(double rho, const std::vector& x, const std::vector>& d_BH, const std::vector>& lambda); + std::vector> dBdrho_func(double rho, double T, const std::vector& x, + const std::vector>& lambda); std::vector> a1ij_func(double rho, double T, const std::vector& x); @@ -152,6 +172,13 @@ class MieKinGas : public Spherical { std::vector> da2ij_div_chi_drho_func(double rho, const std::vector& x, double K_HS, const std::vector>& d_BH, const std::vector>& x0); + std::vector> da2ij_div_chi_drho_func(double rho, double T, const std::vector& x){ + std::vector> d_BH = get_BH_diameters(T); + double zeta_x = zeta_x_func(rho, x, d_BH); + double K_HS = K_HS_func(zeta_x); + std::vector> x0 = get_x0(d_BH); + return da2ij_div_chi_drho_func(rho, x, K_HS, d_BH, x0); + } std::vector> rdf_chi_func(double rho, const std::vector& x, const std::vector>& d_BH); @@ -159,6 +186,11 @@ class MieKinGas : public Spherical { const std::vector>& d_BH); std::vector> gamma_corr(double zeta_x, double T); + std::vector> gamma_corr(double rho, double T, const std::vector& x){ + std::vector> d_BH = get_BH_diameters(T); + double zeta_x = zeta_x_func(rho, x, d_BH); + return gamma_corr(zeta_x, T); + } inline double K_HS_func(double zeta_x){ return pow(1 - zeta_x, 4) / (1 + 4 * zeta_x + 4 * pow(zeta_x, 2) - 4 * pow(zeta_x, 3) + pow(zeta_x, 4)); @@ -191,16 +223,26 @@ class MieKinGas : public Spherical { namespace mie_rdf_constants{ // Gauss Legendre points for computing barker henderson diamenter (see: ThermoPack, SAFT-VR-Mie docs) -constexpr double gl_x[10] = {-0.973906528517171720078, -0.8650633666889845107321, - -0.6794095682990244062343, -0.4333953941292471907993, - -0.1488743389816312108848, 0.1488743389816312108848, - 0.4333953941292471907993, 0.6794095682990244062343, - 0.8650633666889845107321, 0.973906528517171720078}; -constexpr double gl_w[10] = {0.0666713443086881375936, 0.149451349150580593146, - 0.219086362515982043996, 0.2692667193099963550912, - 0.2955242247147528701739, 0.295524224714752870174, - 0.269266719309996355091, 0.2190863625159820439955, - 0.1494513491505805931458, 0.0666713443086881375936}; +constexpr double gl_x[20] = {-0.9931285991850949, -0.9639719272779139, + -0.912234428251326, -0.8391169718222187, + -0.7463319064601508, -0.636053680726515, + -0.5108670019508271, -0.37370608871541955, + -0.22778585114164507, -0.07652652113349737, + 0.07652652113349737, 0.22778585114164507, + 0.37370608871541955, 0.5108670019508271, + 0.636053680726515, 0.7463319064601508, + 0.8391169718222187, 0.912234428251326, + 0.9639719272779139, 0.9931285991850949}; +constexpr double gl_w[20] = {0.017614007139152742, 0.040601429800386134, + 0.06267204833410799, 0.08327674157670514, + 0.1019301198172403, 0.11819453196151845, + 0.13168863844917675, 0.14209610931838218, + 0.149172986472604, 0.15275338713072611, + 0.15275338713072611, 0.149172986472604, + 0.14209610931838218, 0.13168863844917675, + 0.11819453196151845, 0.1019301198172403, + 0.08327674157670514, 0.06267204833410799, + 0.040601429800386134, 0.017614007139152742}; constexpr double C_coeff_matr[4][4] // See Eq. A18 of J. Chem. Phys. 139, 154504 (2013); https://doi.org/10.1063/1.4819786 { diff --git a/cpp/PseudoHardSphere.h b/cpp/PseudoHardSphere.h index e610ada..089fb49 100644 --- a/cpp/PseudoHardSphere.h +++ b/cpp/PseudoHardSphere.h @@ -59,5 +59,5 @@ class PseudoHardSphere : public Spherical { return rdf; } - std::vector> get_contact_diameters(double rho, double T, const std::vector& x) override {return sigma;} + std::vector> get_collision_diameters(double rho, double T, const std::vector& x) override {return sigma;} }; \ No newline at end of file diff --git a/cpp/Spherical.cpp b/cpp/Spherical.cpp index 5f342db..fcb3ba0 100644 --- a/cpp/Spherical.cpp +++ b/cpp/Spherical.cpp @@ -14,7 +14,7 @@ Spherical::Spherical(std::vector mole_weights, } -double Spherical::omega(const int& i, const int& j, const int& l, const int& r, const double& T){ +double Spherical::omega(int i, int j, int l, int r, double T){ OmegaPoint point{i, j, l, r, T}, sympoint{j, i, l, r, T}; const std::map::iterator pos = omega_map.find(point); if (pos == omega_map.end()){ @@ -29,7 +29,31 @@ double Spherical::omega(const int& i, const int& j, const int& l, const int& r, return pos->second; } -double Spherical::w_integral(const int& i, const int& j, const double& T, const int& l, const int& r){ +double Spherical::omega_tester(int i, int j, int l, int r, double T, IntegrationParam& param){ + double w = w_integral_tester(i, j, T, l, r, param); + if (i == j) return pow(sigma[i][j], 2) * sqrt((PI * BOLTZMANN * T) / m[i]) * w; + return 0.5 * pow(sigma[i][j], 2) * sqrt(2 * PI * BOLTZMANN * T / (m0[i][j] * M[i][j] * M[j][i])) * w; +} + +double Spherical::w_integral_tester(int i, int j, double T, int l, int r, IntegrationParam& param){ + double I = integrate2d(param.origin, param.end, + param.dg, param.db, + param.refinement_levels_g, param.refinement_levels_b, + param.subdomain_dblder_limit, + i, j, T, l, r, + w_integrand_export); + + return I; +} + +double Spherical::w_integral(int i, int j, double T, int l, int r){ + /* + Evaulate the dimensionless collision integral + + See: The Kinetic Gas Theory of Mie fluids (V. G. Jervell, Norwegian University of Science and Technology, 2022) + https://ntnuopen.ntnu.no/ntnu-xmlui/handle/11250/3029213 + For details on the integration routine implemented in integrate2d. + */ Point origin{1e-7, 1e-7}; Point end{8, 5}; double dg{0.5}, db{0.03125}; @@ -47,17 +71,90 @@ double Spherical::w_integral(const int& i, const int& j, const double& T, const return I; } -double Spherical::w_integrand(const int& i, const int& j, const double& T, - const double& g, const double& b, - const int& l, const int& r){ // Using b = b / sigma to better scale the axes. Multiply the final integral by sigma. +double Spherical::w_integrand(int i, int j, double T, + double g, double b, + int l, int r){ // Using b = b / sigma to better scale the axes. Multiply the final integral by sigma. const double chi_val = chi(i, j, T, g, b * sigma[i][j]); return 2 * exp(- pow(g, 2)) * pow(g, 2.0 * r + 3.0) * (1 - pow(cos(chi_val), l)) * b; }; +std::vector> Spherical::get_b_max(double T, std::vector>& ierr){ + // bmax[i][j] in units of sigma[i][j] + // The maximum value of the impact parameter at which deflection angle (chi) is positive + const double g_avg = sqrt(4. / PI); // Average dimensionless relative velocity + std::vector> bmax(Ncomps, std::vector(Ncomps, 0)); + double b, db, chi_val; + for (int i = 0; i < Ncomps; i++){ + for (int j = i; j < Ncomps; j++){ + ierr[i][j] = ierr[j][i] = 0; + b = 0.5; + db = 0.1; + chi_val = chi(i, j, T, g_avg, b * sigma[i][j]); + while (abs(db) > 1e-5){ + if ((chi_val < 0 && db > 0) || (chi_val > 0 && db < 0)){ + db *= -0.5; + } + b += db; + chi_val = chi(i, j, T, g_avg, b * sigma[i][j]); + if (b > 10) { // Unable to converge, this failure should be handled in get_b_max(double T) + ierr[i][j] = 1; + break; + } + } + bmax[j][i] = bmax[i][j] = b; + } + } + return bmax; + } + +std::vector> Spherical::get_b_max(double T){ + std::vector> ierr(Ncomps, std::vector(Ncomps, 1)); + std::vector> b_max = get_b_max(T, ierr); + for (size_t i = 0; i < Ncomps; i++){ + for (size_t j = i; j < Ncomps; j++){ + if (ierr[i][j]) { + std::cout << "Could not compute upper integration limit for collision diameter (" << i << ", " << j + << ") using sigma as fallback value." << std::endl; + b_max[j][i] = b_max[i][j] = sigma[i][j]; + } + } + } + return b_max; +} + +std::vector> Spherical::get_collision_diameters(double rho, double T, const std::vector& x){ + // Evaluate the integral of Eq. (40) in RET for Mie fluids (doi: 10.1063/5.0149865) + // Note: For models with is_idealgas==true, zeros are returned, as there is no excluded volume at infinite dilution. + // Weights and nodes for 6-point Gauss-Legendre quadrature + constexpr int n_gl_points = 6; + constexpr double gl_points[n_gl_points] = {-0.93246951, -0.66120939, -0.23861919, 0.23861919, 0.66120939, 0.93246951}; + constexpr double gl_weights[n_gl_points] = {0.17132449, 0.36076157, 0.46791393, 0.46791393, 0.36076157, 0.17132449}; + + const double g_avg = sqrt(1. / PI); // Average dimensionless relative speed of collision + std::vector> avg_R(Ncomps, std::vector(Ncomps, 0.)); + if (is_idealgas) { + return avg_R; + } + std::vector> bmax = get_b_max(T); + double point, weight; + + for (int i = 0; i < Ncomps; i++){ + for (int j = i; j < Ncomps; j++){ + for (int p = 0; p < n_gl_points; p++){ + point = gl_points[p] * (bmax[i][j] / 2.) + bmax[i][j] / 2.; + weight = gl_weights[p] * bmax[i][j] / 2.; + avg_R[i][j] += get_R(i, j, T, g_avg, point * sigma[i][j]) * weight / bmax[i][j]; + } + avg_R[j][i] = avg_R[i][j]; + } + } + return avg_R; +} + /* Contains functions for computing the collision integrals Common variables are: - ij : collision type (1 for like molecules of type 1, 2 for type 2, and 12 or 21 for collisions of unlike molecules + i, j : Species indices T : Temperature [K] l, r : Collision integral indices (in omega and w_ functions) g : Dimentionless relative velocity of molecules @@ -72,7 +169,7 @@ Common variables are: #include "KineticGas.h" #include "Integration/Integration.h" -#pragma region // Helper funcions for computing dimentionless collision integrals +// Helper funcions for computing dimentionless collision integrals double Spherical::theta(int i, int j, const double T, const double g, const double b){ #ifdef DEBUG diff --git a/cpp/Spherical.h b/cpp/Spherical.h index 548ed77..55a29e9 100644 --- a/cpp/Spherical.h +++ b/cpp/Spherical.h @@ -11,39 +11,74 @@ Contains: The abstract class 'Spherical', inheriting from 'KineticGas'. This cla #include "Integration/Integration.h" #include "global_params.h" +class IntegrationParam; + class Spherical : public KineticGas { public: - // In the general case, sigma is a scaling parameter - // On the order of the molecular size. Its specific physical meaning - // is different for different potential models. It may also be without physical meaning - std::vector> sigma; Spherical(std::vector mole_weights, std::vector> sigmaij, bool is_idealgas); virtual ~Spherical(){}; - // Potential models, these must be explicitly overridden in derived classes (in addition to the get_rdf function of KineticGas) + // Potential models, these must be overridden in derived classes, in addition to model_rdf and get_contact_diameters. virtual double potential(int i, int j, double r) = 0; virtual double potential_derivative_r(int i, int j, double r) = 0; virtual double potential_dblderivative_rr(int i, int j, double r) = 0; - double omega(const int& i, const int& j, const int& l, const int& r, const double& T) override; - double w_integral(const int& i, const int& j, const double& T, const int& l, const int& r); // Dimentionless collision integral for spherical potentials - double w_integrand(const int& i, const int& j, const double& T, - const double& g, const double& b, - const int& l, const int& r); - std::function w_integrand_export; // Will bind w_integrand to this function such that it can be passed to the external integration module + double omega(int i, int j, int l, int r, double T) override; + double omega_tester(int i, int j, int l, int r, double T, IntegrationParam& param); + double w_integral_tester(int i, int j, double T, int l, int r, IntegrationParam& param); + std::vector> get_collision_diameters(double rho, double T, const std::vector& x) override; + std::vector> model_rdf(double rho, double T, const std::vector& mole_fracs) override = 0; + + // ------------------------------------------------------------------------------------------------------------------- // + // -------------------------- Spherical Internals are below here ----------------------------------------------------- // + // ---------------------- End users should not need to care about anything below -------------------------------------- // + // ------------------------------------------------------------------------------------------------------------------- // + protected: + // In the general case, sigma is a scaling parameter + // On the order of the molecular size. Its specific physical meaning + // is different for different potential models. It may also be without physical meaning + std::vector> sigma; + + // This method is responsible for calling get_b_max(T, ierr), and handling eventual failures. + // Most inheritance should only require overriding the failure handling, not the computation in get_b_max(T, ierr) + virtual std::vector> get_b_max(double T); + std::vector> get_b_max(double T, std::vector>& ierr); - // Helper functions for computing dimentionless collision integrals double theta(int i, int j, const double T, const double g, const double b); + double chi(int i, int j, double T, double g, double b); + double get_R(int i, int j, double T, double g, double b); + + private: + // Helper functions for computing dimentionless collision integrals double theta_lim(int i, int j, const double T, const double g); double theta_integral(int i, int j, const double T, const double R, const double g, const double b); - virtual double theta_integrand(int i, int j, double T, double r, double g, double b); - virtual double transformed_theta_integrand(int i, int j, double T, double u, double R, double g, double b); - virtual double theta_integrand_dblderivative(int i, int j, double T, double r, double g, double b); - double get_R(int i, int j, double T, double g, double b); - virtual double get_R_rootfunc(int i, int j, double T, double g, double b, double& r); - virtual double get_R_rootfunc_derivative(int i, int j, double T, double g, double b, double& r); - double chi(int i, int j, double T, double g, double b); + double theta_integrand(int i, int j, double T, double r, double g, double b); + double transformed_theta_integrand(int i, int j, double T, double u, double R, double g, double b); + double theta_integrand_dblderivative(int i, int j, double T, double r, double g, double b); + double get_R_rootfunc(int i, int j, double T, double g, double b, double& r); + double get_R_rootfunc_derivative(int i, int j, double T, double g, double b, double& r); + + double w_integral(int i, int j, double T, int l, int r); // Dimentionless collision integral for spherical potentials + double w_integrand(int i, int j, double T, double g, double b, int l, int r); + std::function w_integrand_export; // Will bind w_integrand to this function such that it can be passed to the external integration module + +}; + +class IntegrationParam{ + public: + Point origin{1e-7, 1e-7}; + Point end{8, 5}; + double dg{0.5}, db{0.03125}; + int refinement_levels_g{4}; + int refinement_levels_b{16}; + double subdomain_dblder_limit{1e-5}; + + void set_dg(double dg_){dg = dg_;} + void set_db(double db_){db = db_;} + void set_rlg(int rlg){refinement_levels_g = rlg;} + void set_rlb(int rlb){refinement_levels_b = rlb;} + void set_dblder_lim(double lim){subdomain_dblder_limit = lim;} }; \ No newline at end of file diff --git a/cpp/bindings.cpp b/cpp/bindings.cpp index 9dc5eda..4bea8c3 100644 --- a/cpp/bindings.cpp +++ b/cpp/bindings.cpp @@ -9,32 +9,37 @@ #include namespace py = pybind11; +using vector1d = std::vector; +using vector2d = std::vector; #define KineticGas_bindings(Model) \ - .def("A", &Model::A) \ - .def("A_prime", &Model::A_prime) \ - \ - .def("H_ij", &Model::H_ij) \ - .def("H_i", &Model::H_i) \ .def("get_conductivity_vector", &Model::get_conductivity_vector) \ .def("get_diffusion_vector", &Model::get_diffusion_vector) \ .def("get_diffusion_matrix", &Model::get_diffusion_matrix) \ .def("get_conductivity_matrix", &Model::get_conductivity_matrix) \ .def("get_viscosity_matrix", &Model::get_viscosity_matrix)\ .def("get_viscosity_vector", &Model::get_viscosity_vector)\ - \ - .def("B_prime", &Model::B_prime) \ - .def("B_dblprime", &Model::B_dblprime) \ - \ - .def("L_ij", &Model::L_ij) \ - .def("L_i", &Model::L_i) \ - \ + .def("get_collision_diameters", &Model::get_collision_diameters) \ .def("get_rdf", &Model::get_rdf) \ .def("get_K_factors", &Model::get_K_factors) \ - .def("get_K_prime_factors", &Model::get_K_prime_factors)\ - .def("get_contact_diameters", &Model::get_contact_diameters) \ - \ - .def_readwrite("omega_map", &Model::omega_map) + .def("get_K_prime_factors", &Model::get_K_prime_factors) + /* + .def_readwrite("omega_map", &Model::omega_map) + .def("A", &Model::A) \ + .def("A_prime", &Model::A_prime) \ + \ + .def("H_ij", &Model::H_ij) \ + .def("H_i", &Model::H_i) \ + \ + .def("B_prime", &Model::B_prime) \ + .def("B_dblprime", &Model::B_dblprime) \ + \ + .def("L_ij", &Model::L_ij) \ + .def("L_i", &Model::L_i) \ + \ + + */ + #define Spherical_potential_bindings(Model) \ .def("potential", &Model::potential) \ @@ -42,32 +47,23 @@ namespace py = pybind11; .def("potential_dblderivative_rr", &Model::potential_dblderivative_rr) \ #define Spherical_bindings(Model) \ - .def("chi", &Model::chi) \ - .def("get_R", &Model::get_R) \ - .def("omega", &Model::omega) \ - \ - .def("get_R_rootfunc", &Model::get_R_rootfunc) \ - .def("get_R_rootfunc_derivative", &Model::get_R_rootfunc_derivative) \ - \ - .def("theta", &Model::theta) \ - .def("theta_lim", &Model::theta_lim) \ - .def("theta_integral", &Model::theta_integral) \ - .def("theta_integrand", &Model::theta_integrand) \ - .def("theta_integrand_dblderivative", &Model::theta_integrand_dblderivative) \ - \ - .def("w_integrand", &Model::w_integrand) \ - .def("w_integral", &Model::w_integral) - -#define Mie_rdf_bindings(Model) \ - .def("model_rdf", py::overload_cast&>(&Model::model_rdf))\ - .def("gHS", &Model::rdf_HS)\ - .def("g1", py::overload_cast&, \ - const std::vector>&>(&Model::rdf_g1_func))\ - .def("g2", py::overload_cast&,\ - const std::vector>&, \ - const std::vector>&>(&Model::rdf_g2_func)) + .def("omega", &Model::omega) + /* + .def("chi", &Model::chi) \ + .def("get_R", &Model::get_R) \ + \ + .def("get_R_rootfunc", &Model::get_R_rootfunc) \ + .def("get_R_rootfunc_derivative", &Model::get_R_rootfunc_derivative) \ + \ + .def("theta", &Model::theta) \ + .def("theta_lim", &Model::theta_lim) \ + .def("theta_integral", &Model::theta_integral) \ + .def("theta_integrand", &Model::theta_integrand) \ + .def("theta_integrand_dblderivative", &Model::theta_integrand_dblderivative) \ + \ + .def("w_integrand", &Model::w_integrand) \ + .def("w_integral", &Model::w_integral) + */ #ifndef DEBUG @@ -79,11 +75,6 @@ PYBIND11_MODULE(KineticGas_d, handle){ handle.def("ipow", &ipow); handle.def("factorial_tests", &factorial_tests); - // handle.def("zeta_x", &mie_rdf::zeta_x_func); - // handle.def("dzeta_x_drho", &mie_rdf::dzetax_drho_func); - // handle.def("zeta_eff", &mie_rdf::zeta_eff_func); - // handle.def("dzeta_eff_drho", &mie_rdf::dzeta_eff_drho_func); - py::class_(handle, "Product") .def(py::init()) .def(py::init()) @@ -108,38 +99,49 @@ PYBIND11_MODULE(KineticGas_d, handle){ }); py::class_(handle, "cpp_MieKinGas") - .def(py::init< - std::vector, - std::vector>, - std::vector>, - std::vector>, - std::vector>, - bool + .def(py::init() ) KineticGas_bindings(MieKinGas) Spherical_potential_bindings(MieKinGas) Spherical_bindings(MieKinGas) - Mie_rdf_bindings(MieKinGas) - .def("g1", py::overload_cast&>(&MieKinGas::rdf_g1_func)) - .def("g2", py::overload_cast&>(&MieKinGas::rdf_g2_func)) - .def("a1s", py::overload_cast&, - const std::vector>&>(&MieKinGas::a_1s_func)) - .def("da1s_drho", py::overload_cast&, - const std::vector>&>(&MieKinGas::da1s_drho_func)) - .def("get_dBH", &MieKinGas::get_BH_diameters) - .def("a1ij", &MieKinGas::a1ij_func) - .def("da1ij_drho", py::overload_cast&>(&MieKinGas::da1ij_drho_func)) - .def("a2ij", py::overload_cast&>(&MieKinGas::a2ij_func)) - .def("da2ij_drho", py::overload_cast&>(&MieKinGas::da2ij_drho_func)) - .def("get_BH_diameters", &MieKinGas::get_BH_diameters); + .def("get_BH_diameters", &MieKinGas::get_BH_diameters) + .def("get_vdw_alpha", &MieKinGas::get_vdw_alpha) + .def("saft_rdf", &MieKinGas::saft_rdf) + .def("rdf_g0", py::overload_cast(&MieKinGas::rdf_HS)) + .def("rdf_g1", py::overload_cast(&MieKinGas::rdf_g1_func)) + .def("rdf_g2", py::overload_cast(&MieKinGas::rdf_g2_func)) + .def("a1", &MieKinGas::a1ij_func) + .def("da1drho", py::overload_cast(&MieKinGas::da1ij_drho_func)) + .def("a1s", py::overload_cast(&MieKinGas::a_1s_func)) + .def("da1s_drho", py::overload_cast(&MieKinGas::da1s_drho_func)) + .def("B_func", py::overload_cast(&MieKinGas::B_func)) + .def("dBdrho_func", py::overload_cast(&MieKinGas::dBdrho_func)) + ; + // .def("da1_drho", py::overload_cast(&MieKinGas::da1ij_drho_func)) + // .def("da1s_drho", py::overload_cast(&MieKinGas::da1s_drho_func)) + // .def("zeta_eff", &MieKinGas::zeta_eff_func) + // .def("dzeta_eff_drho", &MieKinGas::dzeta_eff_drho_func) + // .def("zeta_x", &MieKinGas::zeta_x_func) + // .def("a1s", py::overload_cast(&MieKinGas::a_1s_func)) + // .def("B_func", py::overload_cast(&MieKinGas::B_func)) + // .def("get_dBH", &MieKinGas::get_BH_diameters) + // .def("a1ij", &MieKinGas::a1ij_func) + // .def("a2ij", py::overload_cast(&MieKinGas::a2ij_func)) + // .def("da2ij_drho", py::overload_cast(&MieKinGas::da2ij_drho_func)) + // .def("da2_div_chi_drho", py::overload_cast(&MieKinGas::da2ij_div_chi_drho_func)) + // .def("gamma_corr", py::overload_cast(&MieKinGas::gamma_corr)) + // ; py::class_(handle, "cpp_HardSphere") .def(py::init< - std::vector, - std::vector>, + vector1d, + vector2d, bool >() ) @@ -151,8 +153,8 @@ PYBIND11_MODULE(KineticGas_d, handle){ py::class_(handle, "cpp_PseudoHardSphere") .def(py::init< - std::vector, - std::vector>, + vector1d, + vector2d, bool >() ) diff --git a/cpp/build.sh b/cpp/build.sh index 8aa2e24..41248e6 100644 --- a/cpp/build.sh +++ b/cpp/build.sh @@ -2,7 +2,7 @@ # Script to build the KineticGas module and run tests. Can be used with options set -e -py_version=${1-"-DPYBIND11_PYTHON_VERSION=3"} +py_version=${1-"-DPYBIND11_PYTHON_VERSION=3 -DPYTHON_EXECUTABLE=$(which python)"} export CC=/opt/homebrew/bin/gcc-13 # Modify as needed, commonly set to /usr/bin/gcc export CXX=/opt/homebrew/bin/g++-13 # Modify as needed, commonly set to /usr/bin/g++ diff --git a/docs/vCurrent/HardSphere_methods.md b/docs/vCurrent/HardSphere_methods.md index 0371ea3..bdf6b21 100644 --- a/docs/vCurrent/HardSphere_methods.md +++ b/docs/vCurrent/HardSphere_methods.md @@ -6,7 +6,7 @@ permalink: /vcurrent/HardSphere_methods.html ---