diff --git a/include/teqp/derivs.hpp b/include/teqp/derivs.hpp
index 0233b2ba78e87cf631321a12c79825b346450a7b..fa1a81ed091024bb76ff41fd1c2c5b2b0bd9a8e5 100644
--- a/include/teqp/derivs.hpp
+++ b/include/teqp/derivs.hpp
@@ -67,9 +67,28 @@ enum class ADBackends { autodiff, multicomplex, complex_step };
template<typename Model, typename Scalar = double, typename VectorType = Eigen::ArrayXd>
struct TDXDerivatives {
+ template<ADBackends be = ADBackends::autodiff>
static auto get_Ar10(const Model& model, const Scalar &T, const Scalar &rho, const VectorType& molefrac) {
- double h = 1e-100;
- return -T*model.alphar(std::complex<Scalar>(T, h), rho, molefrac).imag()/h; // Complex step derivative
+ if constexpr (be == ADBackends::complex_step) {
+ double h = 1e-100;
+ return -T * model.alphar(std::complex<Scalar>(T, h), rho, molefrac).imag() / h; // Complex step derivative
+ }
+ else if constexpr (be == ADBackends::multicomplex) {
+ using fcn_t = std::function<mcx::MultiComplex<Scalar>(const mcx::MultiComplex<Scalar>&)>;
+ bool and_val = true;
+ fcn_t f = [&model, &rho, &molefrac](const auto& Trecip_) { return model.alphar(1.0/Trecip_, rho, molefrac); };
+ auto ders = diff_mcx1(f, 1.0/T, 1, and_val);
+ return (1.0/T)*ders[1];
+ }
+ else if constexpr (be == ADBackends::autodiff) {
+ autodiff::dual Trecipdual = 1.0/T;
+ auto f = [&model, &rho, &molefrac](const auto& Trecip_) { return eval(model.alphar(eval(1.0/Trecip_), rho, molefrac)); };
+ auto der = derivative(f, wrt(Trecipdual), at(Trecipdual));
+ return (1.0/T)*der;
+ }
+ else {
+ static_assert("algorithmic differentiation backend is invalid in get_Ar10");
+ }
}
template<ADBackends be = ADBackends::autodiff>
@@ -77,28 +96,85 @@ struct TDXDerivatives {
if constexpr(be == ADBackends::complex_step){
double h = 1e-100;
auto der = model.alphar(T, std::complex<Scalar>(rho, h), molefrac).imag() / h;
- return der*rho;
+ return rho*der;
}
else if constexpr(be == ADBackends::multicomplex){
using fcn_t = std::function<mcx::MultiComplex<Scalar>(const mcx::MultiComplex<Scalar>&)>;
bool and_val = true;
fcn_t f = [&model, &T, &molefrac](const auto& rho_) { return model.alphar(T, rho_, molefrac); };
auto ders = diff_mcx1(f, rho, 1, and_val);
- return ders[1] * rho;
+ return rho*ders[1];
}
else if constexpr(be == ADBackends::autodiff){
autodiff::dual rhodual = rho;
auto f = [&model, &T, &molefrac](const auto& rho_) { return eval(model.alphar(T, rho_, molefrac)); };
auto der = derivative(f, wrt(rhodual), at(rhodual));
- return der * rho;
+ return rho*der;
+ }
+ else {
+ static_assert("algorithmic differentiation backend is invalid in get_Ar01");
}
}
+ template<ADBackends be = ADBackends::autodiff>
static auto get_Ar02(const Model& model, const Scalar& T, const Scalar& rho, const VectorType& molefrac) {
- autodiff::dual2nd rhodual = rho;
- auto f = [&model, &T, &molefrac](const auto& rho_) { return eval(model.alphar(T, rho_, molefrac)); };
- auto ders = derivatives(f, wrt(rhodual), at(rhodual));
- return ders[2]*rho*rho;
+ if constexpr (be == ADBackends::autodiff) {
+ autodiff::dual2nd rhodual = rho;
+ auto f = [&model, &T, &molefrac](const auto& rho_) { return eval(model.alphar(T, rho_, molefrac)); };
+ auto ders = derivatives(f, wrt(rhodual), at(rhodual));
+ return rho*rho*ders[2];
+ }
+ else {
+ static_assert("algorithmic differentiation backend is invalid in get_Ar02");
+ }
+ }
+
+ template<ADBackends be = ADBackends::autodiff>
+ static auto get_Ar20(const Model& model, const Scalar& T, const Scalar& rho, const VectorType& molefrac) {
+ if constexpr (be == ADBackends::autodiff) {
+ autodiff::dual2nd Trecipdual = 1/T;
+ auto f = [&model, &rho, &molefrac](const auto& Trecip) { return eval(model.alphar(eval(1/Trecip), rho, molefrac)); };
+ auto ders = derivatives(f, wrt(Trecipdual), at(Trecipdual));
+ return (1/T)*(1/T)*ders[2];
+ }
+ else {
+ static_assert("algorithmic differentiation backend is invalid in get_Ar20");
+ }
+ }
+
+ template<ADBackends be = ADBackends::autodiff>
+ static auto get_Ar11(const Model& model, const Scalar& T, const Scalar& rho, const VectorType& molefrac) {
+
+ if constexpr (be == ADBackends::multicomplex) {
+ using fcn_t = std::function< mcx::MultiComplex<double>(const std::valarray<mcx::MultiComplex<double>>&)>;
+ const fcn_t func = [&model, &molefrac](const auto& zs) {
+ auto rhomolar = zs[0], Trecip = zs[1];
+ return model.alphar(1.0/Trecip, rhomolar, molefrac);
+ };
+ std::vector<double> xs = { rho, 1.0/T};
+ std::vector<int> order = { 1, 1 };
+ auto der = mcx::diff_mcxN(func, xs, order);
+ return (1.0/T)*rho*der;
+ }
+ else if constexpr (be == ADBackends::autodiff) {
+ //static_assert("bug in autodiff, can't use autodiff for cross derivative");
+ autodiff::dual2nd rhodual = rho, Trecipdual=1/T;
+ auto f = [&model, &molefrac](const auto& Trecip, const auto&rho_) { return eval(model.alphar(eval(1/Trecip), rho_, molefrac)); };
+ //auto der = derivative(f, wrt(Trecipdual, rhodual), at(Trecipdual, rhodual)); // d^2alphar/drhod(1/T) // This should work, but gives instead 1,0 derivative
+ auto [u01, u10, u11] = derivatives(f, wrt(Trecipdual, rhodual), at(Trecipdual, rhodual)); // d^2alphar/drhod(1/T)
+ return (1.0/T)*rho*u11;
+ }
+ else {
+ static_assert("algorithmic differentiation backend is invalid in get_Ar11");
+ }
+ }
+
+ template<ADBackends be = ADBackends::autodiff>
+ static auto get_neff(const Model& model, const Scalar& T, const Scalar& rho, const VectorType& molefrac) {
+ auto Ar01 = get_Ar01<be>(model, T, rho, molefrac);
+ auto Ar11 = get_Ar11<be>(model, T, rho, molefrac);
+ auto Ar20 = get_Ar20(model, T, rho, molefrac);
+ return -3*(Ar01-Ar11)/Ar20;
}
};
diff --git a/include/teqp/models/eos.hpp b/include/teqp/models/eos.hpp
index e0dff36b423830335cf0d38cdba462e8e2e21f3a..c9e77c2c9e05886990254b8c0f9846019a64d9d0 100644
--- a/include/teqp/models/eos.hpp
+++ b/include/teqp/models/eos.hpp
@@ -55,7 +55,7 @@ protected:
a_ = a_ + molefracs[i] * molefracs[j] * aij;
}
}
- return a_;
+ return forceeval(a_);
}
template<typename CompType>
@@ -64,7 +64,7 @@ protected:
for (auto i = 0; i < molefracs.size(); ++i) {
b_ = b_ + molefracs[i] * bi[i];
}
- return b_;
+ return forceeval(b_);
}
public:
diff --git a/include/teqp/models/pcsaft.hpp b/include/teqp/models/pcsaft.hpp
index 4e9f7e2a17568e06294b5a6c3d11d3cf61466f88..b1152f894607fcc4aed5819fed7d6cb5a6cc011f 100644
--- a/include/teqp/models/pcsaft.hpp
+++ b/include/teqp/models/pcsaft.hpp
@@ -261,6 +261,6 @@ public:
}
auto alphar_hc = mbar * get_alphar_hs(zeta) - sumproduct(mole_fractions, mminus1, lngii_hs); // Eq. A.4
auto alphar_disp = -2 * MY_PI * rho_A3 * I1 * c.m2_epsilon_sigma3_bar - MY_PI * rho_A3 * mbar * C1(eta, mbar) * I2 * c.m2_epsilon2_sigma3_bar;
- return alphar_hc + alphar_disp;
+ return forceeval(alphar_hc + alphar_disp);
}
};
\ No newline at end of file
diff --git a/interface/pybind11_wrapper.cpp b/interface/pybind11_wrapper.cpp
index ace92566e58dc9fca20dd23ab40f6c8279a53fc9..b42694b238e1b22bc5ed4a5ca8197e4f31af6afc 100644
--- a/interface/pybind11_wrapper.cpp
+++ b/interface/pybind11_wrapper.cpp
@@ -13,9 +13,12 @@ template<typename Model>
void add_TDx_derivatives(py::module& m) {
using id = TDXDerivatives<Model, double, Eigen::Array<double, Eigen::Dynamic, 1> >;
//m.def("get_Ar00", &id::get_Ar00, py::arg("model"), py::arg("T"), py::arg("rho"), py::arg("molefrac"));
- m.def("get_Ar10", &id::get_Ar10, py::arg("model"), py::arg("T"), py::arg("rho"), py::arg("molefrac"));
+ m.def("get_Ar10", &id::get_Ar10<ADBackends::autodiff>, py::arg("model"), py::arg("T"), py::arg("rho"), py::arg("molefrac"));
m.def("get_Ar01", &id::get_Ar01<ADBackends::autodiff>, py::arg("model"), py::arg("T"), py::arg("rho"), py::arg("molefrac"));
- m.def("get_Ar02", &id::get_Ar02, py::arg("model"), py::arg("T"), py::arg("rho"), py::arg("molefrac"));
+ m.def("get_Ar11", &id::get_Ar11<ADBackends::autodiff>, py::arg("model"), py::arg("T"), py::arg("rho"), py::arg("molefrac"));
+ m.def("get_Ar02", &id::get_Ar02<ADBackends::autodiff>, py::arg("model"), py::arg("T"), py::arg("rho"), py::arg("molefrac"));
+ m.def("get_Ar20", &id::get_Ar20<ADBackends::autodiff>, py::arg("model"), py::arg("T"), py::arg("rho"), py::arg("molefrac"));
+ m.def("get_neff", &id::get_neff<ADBackends::autodiff>, py::arg("model"), py::arg("T"), py::arg("rho"), py::arg("molefrac"));
}
template<typename Model>
diff --git a/src/main.cpp b/src/main.cpp
index 25e513040f61ceea037962fc16cd2c02857be3c1..c4e218f067ef40737f317d05b354a00e510a1512 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -8,74 +8,67 @@
#include "teqp/derivs.hpp"
#include "teqp/models/eos.hpp"
#include "teqp/models/pcsaft.hpp"
-//
-//void test_vdwMix() {
-// // Argon + Xenon
-// std::valarray<double> Tc_K = { 150.687, 289.733 };
-// std::valarray<double> pc_Pa = { 4863000.0, 5842000.0 };
-// vdWEOS<double> vdW(Tc_K, pc_Pa);
-//
-// double T = 298.15;
-// auto rho = 3.0;
-// auto R = get_R_gas<double>();
-// auto rhotot = rho;
-//
-// const auto rhovec = (Eigen::ArrayXd(2) << rho/2, rho/2).finished();
-//
-// auto fPsir = [&vdW](const auto& T, const auto& rhovec) {
-// using container = decltype(rhovec);
-// auto rhotot_ = std::accumulate(std::begin(rhovec), std::end(rhovec), (decltype(rhovec[0]))0.0);
-// auto molefrac = rhovec/rhotot_;
-// return vdW.alphar(T, rhotot_, molefrac)*vdW.R*T*rhotot_;
-// };
-// auto Psir = fPsir(T, rhovec);
-// auto dPsirdrho0 = rhovec[0]*derivrhoi(fPsir, T, rhovec, 0);
-// auto dPsirdrho1 = rhovec[1]*derivrhoi(fPsir, T, rhovec, 1);
-// auto pfromderiv = rho*R*T - Psir + dPsirdrho0 + dPsirdrho1;
-// using id = IsochoricDerivatives<decltype(vdW)>;
-// auto sr = id::get_splus(vdW, T, rhovec);
-//
-// auto t2 = std::chrono::steady_clock::now();
-// auto pfromderiv3 = rhotot*R*T + id::get_pr(vdW, T, rhovec);
-// auto t3 = std::chrono::steady_clock::now();
-// std::cout << std::chrono::duration<double>(t3 - t2).count() << " from isochoric (mix) " << std::endl;
-//
-//}
-template<typename T1, typename T2, typename T3>
-void f(const T1& v1, const T2& v2, const T3& v3) {
- using t = decltype(forceeval(v1* v2* v3[0]));
- std::cout << "Hi";
+void test_vdwMix() {
+ // Argon + Xenon
+ std::valarray<double> Tc_K = { 150.687, 289.733 };
+ std::valarray<double> pc_Pa = { 4863000.0, 5842000.0 };
+ vdWEOS<double> vdW(Tc_K, pc_Pa);
+ double T = 298.15;
+ auto rho = 3.0;
+ auto R = get_R_gas<double>();
+ auto rhotot = rho;
+
+ const auto rhovec = (Eigen::ArrayXd(2) << rho/2, rho/2).finished();
+
+ auto fPsir = [&vdW](const auto& T, const auto& rhovec) {
+ using container = decltype(rhovec);
+ auto rhotot_ = std::accumulate(std::begin(rhovec), std::end(rhovec), (decltype(rhovec[0]))0.0);
+ auto molefrac = rhovec/rhotot_;
+ return vdW.alphar(T, rhotot_, molefrac)*vdW.R*T*rhotot_;
+ };
+ auto Psir = fPsir(T, rhovec);
+ auto dPsirdrho0 = rhovec[0]*derivrhoi(fPsir, T, rhovec, 0);
+ auto dPsirdrho1 = rhovec[1]*derivrhoi(fPsir, T, rhovec, 1);
+ auto pfromderiv = rho*R*T - Psir + dPsirdrho0 + dPsirdrho1;
+ using id = IsochoricDerivatives<decltype(vdW)>;
+ auto sr = id::get_splus(vdW, T, rhovec);
+
+ auto t2 = std::chrono::steady_clock::now();
+ auto pfromderiv3 = rhotot*R*T + id::get_pr(vdW, T, rhovec);
+ auto t3 = std::chrono::steady_clock::now();
+ std::cout << std::chrono::duration<double>(t3 - t2).count() << " from isochoric (mix) " << std::endl;
}
int main(){
- //test_vdwMix();
-
- //autodiff::dual1st x;
- //double y;
- //Eigen::VectorX<double> z;
- //f(x, y, z);
- ////using t = decltype(forceeval(x * y * z[0]));
- ////t f = "";
+ test_vdwMix();
std::vector<std::string> names = { "Methane", "Ethane" };
PCSAFTMixture mix(names);
mix.print_info();
using id = IsochoricDerivatives<decltype(mix)>;
using vd = VirialDerivatives<decltype(mix)>;
+ using tdx = TDXDerivatives<decltype(mix)>;
double T = 300;
const auto rhovec = (Eigen::ArrayXd(2) << 1.0, 2.0).finished();
const Eigen::ArrayXd molefrac = (rhovec/rhovec.sum()).eval();
const double rho = rhovec.sum();
- double a00csd = get_Ar01<ADBackends::complex_step>(mix, T, rho, molefrac);
- double a00cx = get_Ar01<ADBackends::multicomplex>(mix, T, rho, molefrac);
- double a00ad = get_Ar01<ADBackends::autodiff>(mix, T, rho, molefrac);
+
+ double a00csd = tdx::get_Ar01<ADBackends::complex_step>(mix, T, rho, molefrac);
+ double a00cx = tdx::get_Ar01<ADBackends::multicomplex>(mix, T, rho, molefrac);
+ double a00ad = tdx::get_Ar01<ADBackends::autodiff>(mix, T, rho, molefrac);
+ double a11ad = tdx::get_Ar11<ADBackends::autodiff>(mix, T, rho, molefrac);
+
+ for (auto T = 20; T < 1e5; T *= 1.3) {
+ double neff = tdx::get_neff<ADBackends::autodiff>(mix, T, 1e-10, molefrac);
+ std::cout << T << "," << neff << std::endl;
+ }
+ double neff = tdx::get_neff<ADBackends::autodiff>(mix, T, rho, molefrac);
double a00iso = id::get_Ar01(mix, T, rhovec);
double apriso = id::get_pr(mix, T, rhovec);
double B2 = vd::get_B2vir(mix, T, molefrac);
double B12 = vd::get_B12vir(mix, T, molefrac);
-
return EXIT_SUCCESS;
}
\ No newline at end of file
diff --git a/src/multifluid.cpp b/src/multifluid.cpp
index b180c93463f4bf5c520304c9668dd1455795de45..a8abfd0b9131849b63ed782d56d496711237a154 100644
--- a/src/multifluid.cpp
+++ b/src/multifluid.cpp
@@ -62,30 +62,17 @@ void trace_critical_loci(const std::string &coolprop_root, const nlohmann::json
}
}
-int main(){
-
- std::string coolprop_root = "C:/Users/ihb/Code/CoolProp";
- coolprop_root = "../mycp";
- auto BIPcollection = nlohmann::json::parse(
- std::ifstream(coolprop_root + "/dev/mixtures/mixture_binary_pairs.json")
- );
-
- //// Critical curves
- //{
- // Timer t(1);
- // trace_critical_loci(coolprop_root, BIPcollection);
- //}
-
-{
+template<typename J>
+void time_calls(const std::string &coolprop_root, const J &BIPcollection) {
auto model = build_multifluid_model({ "methane", "ethane" }, coolprop_root, BIPcollection);
- Eigen::ArrayXd rhovec(2); rhovec << 1.0, 2.0 ;
+ Eigen::ArrayXd rhovec(2); rhovec << 1.0, 2.0;
double T = 300;
{
- const auto molefrac = (Eigen::ArrayXd(2) << rhovec[0]/rhovec.sum(), rhovec[1]/rhovec.sum()).finished();
+ const auto molefrac = (Eigen::ArrayXd(2) << rhovec[0] / rhovec.sum(), rhovec[1] / rhovec.sum()).finished();
using vd = VirialDerivatives<decltype(model)>;
auto B12 = vd::get_B12vir(model, T, molefrac);
-
+
using id = IsochoricDerivatives<decltype(model)>;
auto mu = id::get_chempot_autodiff(model, T, rhovec);
@@ -93,11 +80,12 @@ int main(){
double T = 300.0;
constexpr int N = 10000;
volatile double alphar;
- double rrrr = get_Ar01(model, T, rho, molefrac);
- double rrrr2 = get_Ar02(model, T, rho, molefrac);
+ using tdx = TDXDerivatives<decltype(model)>;
+ double rrrr = tdx::get_Ar01(model, T, rho, molefrac);
+ double rrrr2 = tdx::get_Ar02(model, T, rho, molefrac);
{
Timer t(N);
- for (auto i = 0; i < N; ++i){
+ for (auto i = 0; i < N; ++i) {
alphar = model.alphar(T, rho, molefrac);
}
std::cout << alphar << " function call" << std::endl;
@@ -105,28 +93,28 @@ int main(){
{
Timer t(N);
for (auto i = 0; i < N; ++i) {
- alphar = get_Ar01<ADBackends::complex_step>(model, T, rho, molefrac);
+ alphar = tdx::get_Ar01<ADBackends::complex_step>(model, T, rho, molefrac);
}
std::cout << alphar << "; 1st CSD" << std::endl;
}
{
Timer t(N);
for (auto i = 0; i < N; ++i) {
- alphar = get_Ar01<ADBackends::autodiff>(model, T, rho, molefrac);
+ alphar = tdx::get_Ar01<ADBackends::autodiff>(model, T, rho, molefrac);
}
std::cout << alphar << "; 1st autodiff::autodiff" << std::endl;
}
{
Timer t(N);
for (auto i = 0; i < N; ++i) {
- alphar = get_Ar01<ADBackends::multicomplex>(model, T, rho, molefrac);
+ alphar = tdx::get_Ar01<ADBackends::multicomplex>(model, T, rho, molefrac);
}
std::cout << alphar << "; 1st MCX" << std::endl;
}
{
Timer t(N);
for (auto i = 0; i < N; ++i) {
- alphar = get_Ar02(model, T, rho, molefrac);
+ alphar = tdx::get_Ar02(model, T, rho, molefrac);
}
std::cout << alphar << "; 2nd autodiff" << std::endl;
}
@@ -152,14 +140,40 @@ int main(){
std::cout << alphar << "; 5 derivs" << std::endl;
}
}
+}
- const auto molefrac = (Eigen::ArrayXd(2) << 1.0 / 3.0, 2.0 / 3.0 ).finished();
+int main(){
+
+ std::string coolprop_root = "C:/Users/ihb/Code/CoolProp";
+ coolprop_root = "../mycp";
+ auto BIPcollection = nlohmann::json::parse(
+ std::ifstream(coolprop_root + "/dev/mixtures/mixture_binary_pairs.json")
+ );
+
+ //// Critical curves
+ //{
+ // Timer t(1);
+ // trace_critical_loci(coolprop_root, BIPcollection);
+ //}
+
+ //time_calls(coolprop_root, BIPcollection);
+
+{
+ auto model = build_multifluid_model({ "methane", "ethane" }, coolprop_root, BIPcollection);
+ Eigen::ArrayXd rhovec(2); rhovec << 1.0, 2.0;
+ double T = 300;
+ const auto molefrac = rhovec/rhovec.sum();
+ using tdx = TDXDerivatives<decltype(model)>;
+ const auto b = ADBackends::autodiff;
auto rho = rhovec.sum();
auto alphar = model.alphar(T, rho, rhovec);
- auto Ar01 = get_Ar01(model, T, rho, molefrac);
- auto Ar10 = get_Ar10(model, T, rho, molefrac);
- auto Ar02 = get_Ar02(model, T, rho, molefrac);
+ auto Ar01 = tdx::get_Ar01<b>(model, T, rho, molefrac);
+ auto Ar10 = tdx::get_Ar10(model, T, rho, molefrac);
+ auto Ar02 = tdx::get_Ar02(model, T, rho, molefrac);
+ auto Ar11 = tdx::get_Ar11<b>(model, T, rho, molefrac);
+ auto Ar11mcx = tdx::get_Ar11<ADBackends::multicomplex>(model, T, rho, molefrac);
+ auto Ar20 = tdx::get_Ar20(model, T, rho, molefrac);
using id = IsochoricDerivatives<decltype(model)>;
auto splus = id::get_splus(model, T, rhovec);
diff --git a/src/tests/catch_tests.cxx b/src/tests/catch_tests.cxx
index fe6a4e90cec917a8dfc6b597d2c21c89f946ddcf..bd0104c9b838b6c28651506c76d15b055b434e54 100644
--- a/src/tests/catch_tests.cxx
+++ b/src/tests/catch_tests.cxx
@@ -2,6 +2,7 @@
#include "catch/catch.hpp"
#include "teqp/core.hpp"
+#include "teqp/models/pcsaft.hpp"
auto build_vdW_argon() {
double Omega_b = 1.0 / 8, Omega_a = 27.0 / 64;
@@ -83,6 +84,27 @@ TEST_CASE("Check virial coefficients for vdW", "[virial]")
}
}
+TEST_CASE("Check neff", "[virial]")
+{
+ double T = 298.15;
+ double rho = 3.0;
+ const Eigen::Array2d molefrac = { 0.5, 0.5 };
+ auto f = [&T, &rho, &molefrac](const auto& model) {
+ auto neff = TDXDerivatives<decltype(model)>::get_neff(model, T, rho, molefrac);
+ CAPTURE(neff);
+ CHECK(neff > 0);
+ CHECK(neff < 100);
+ };
+ // This quantity is undefined for the van der Waals EOS because Ar20 is always 0
+ //SECTION("vdW") {
+ // f(build_simple());
+ //}
+ SECTION("PCSAFT") {
+ std::vector<std::string> names = { "Methane", "Ethane" };
+ f(PCSAFTMixture(names));
+ }
+}
+
TEST_CASE("Check Hessian of Psir", "[virial]")
{