Skip to content

Commit

Permalink
Merge pull request #24 from thermotools/mie_rdf
Browse files Browse the repository at this point in the history 
Mie rdf
  • Loading branch information
@vegardjervell
vegardjervell authored Apr 16, 2024
2 parents 6eb8338 + 3feca6e commit 8f373b9
Show file tree
Hide file tree
Showing 17 changed files with 529 additions and 307 deletions.
5 changes: 3 additions & 2 deletions cpp/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
@@ -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}")
Expand Down Expand Up @@ -32,4 +33,4 @@ else()
endif()

add_subdirectory(${PYBIND11_ROOT} release)
pybind11_add_module(${TARGET} ${SOURCES})
pybind11_add_module(${TARGET} ${SOURCES})
6 changes: 3 additions & 3 deletions cpp/HardSphere.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -59,6 +59,6 @@ class HardSphere : public KineticGas {
return rdf;
}

std::vector<std::vector<double>> get_contact_diameters(double rho, double T, const std::vector<double>& x) override {return sigma;}
std::vector<std::vector<double>> get_collision_diameters(double rho, double T, const std::vector<double>& x) override {return sigma;}

};
25 changes: 11 additions & 14 deletions cpp/KineticGas.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@ std::vector<double> KineticGas::get_K_factors(double rho, double T, const std::v

std::vector<std::vector<double>> rdf = get_rdf(rho, T, mole_fracs);
std::vector<double> K(Ncomps, 0.);
std::vector<std::vector<double>> cd = get_contact_diameters(rho, T, mole_fracs);
std::vector<std::vector<double>> 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];
Expand All @@ -80,7 +80,7 @@ std::vector<double> KineticGas::get_K_prime_factors(double rho, double T, const

std::vector<std::vector<double>> rdf = get_rdf(rho, T, mole_fracs);
std::vector<double> K_prime(Ncomps, 0.);
std::vector<std::vector<double>> cd = get_contact_diameters(rho, T, mole_fracs);
std::vector<std::vector<double>> 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];
Expand Down Expand Up @@ -293,7 +293,7 @@ std::vector<double> 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++){
Expand All @@ -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;

Expand Down Expand Up @@ -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};
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand All @@ -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};

Expand All @@ -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]};
Expand All @@ -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++){
Expand Down
128 changes: 69 additions & 59 deletions cpp/KineticGas.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<double> m;
std::vector<std::vector<double>> M, m0;
std::map<OmegaPoint, double> omega_map;

KineticGas(std::vector<double> 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<std::vector<double>> get_contact_diameters(double rho, double T, const std::vector<double>& 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<std::vector<double>> model_rdf(double T, double rho, const std::vector<double>& mole_fracs) = 0;
virtual std::vector<std::vector<double>> get_collision_diameters(double rho, double T, const std::vector<double>& x) = 0;

std::vector<double> get_wt_fracs(const std::vector<double> mole_fracs); // Compute weight fractions from mole fractions
// Radial distribution function "at contact". Inheriting classes must implement model_rdf.
std::vector<std::vector<double>> get_rdf(double rho, double T, const std::vector<double>& mole_fracs) {
if (is_idealgas) return std::vector<std::vector<double>>(Ncomps, std::vector<double>(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<std::vector<double>> model_rdf(double rho, double T, const std::vector<double>& 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<double> get_wt_fracs(const std::vector<double> 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<std::vector<double>> get_conductivity_matrix(double rho, double T, const std::vector<double>& x, int N);
std::vector<double> get_diffusion_vector(double rho, double T, const std::vector<double>& x, int N);
Expand All @@ -128,32 +100,70 @@ class KineticGas{
std::vector<double> get_conductivity_vector(double rho, double T, const std::vector<double>& x, int N);
std::vector<std::vector<double>> get_viscosity_matrix(double rho, double T, const std::vector<double>&x, int N);
std::vector<double> get_viscosity_vector(double rho, double T, const std::vector<double>& x, int N);


// Eq. (1.2) of 'multicomponent docs'
std::vector<double> get_K_factors(double rho, double T, const std::vector<double>& mole_fracs);
// Eq. (5.4) of 'multicomponent docs'
std::vector<double> get_K_prime_factors(double rho, double T, const std::vector<double>& mole_fracs);
std::vector<double> get_K_factors(double rho, double T, const std::vector<double>& mole_fracs); // Eq. (1.2) of 'multicomponent docs'
std::vector<double> get_K_prime_factors(double rho, double T, const std::vector<double>& 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<double> m;
std::vector<std::vector<double>> M, m0;
std::map<OmegaPoint, double> 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_<ij>() 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<int>& i_vec, const std::vector<int>& 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<int>& i_vec, const std::vector<int>& j_vec,
const std::vector<int>& l_vec, const std::vector<int>& 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;}
Loading

0 comments on commit 8f373b9

Please sign in to comment.