From c1acbabb18663a254311d3637b63eb3b44f8cb88 Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.bell@nist.gov>
Date: Thu, 1 Apr 2021 09:14:08 -0400
Subject: [PATCH] Add (some) virials to the python interface
---
include/teqp/derivs.hpp | 124 ++++++++++++++++++---------------
interface/pybind11_wrapper.cpp | 12 +++-
src/multifluid.cpp | 11 +--
src/tests/catch_tests.cxx | 17 +++--
4 files changed, 93 insertions(+), 71 deletions(-)
diff --git a/include/teqp/derivs.hpp b/include/teqp/derivs.hpp
index 36ab64e..7b75166 100644
--- a/include/teqp/derivs.hpp
+++ b/include/teqp/derivs.hpp
@@ -62,10 +62,14 @@ typename ContainerType::value_type derivrhoi(const FuncType& f, TType T, const C
return f(T, rhocom).imag() / h;
}
+
+
+
+
template <typename Model, typename TType, typename RhoType, typename ContainerType>
typename ContainerType::value_type get_Ar10(const Model& model, const TType T, const RhoType &rho, const ContainerType& molefrac) {
double h = 1e-100;
- return -T*model.alphar(std::complex<TType>(T, h), rho, molefrac); // Complex step derivative
+ return -T*model.alphar(std::complex<TType>(T, h), rho, molefrac).imag()/h; // Complex step derivative
}
enum class ADBackends { autodiff, multicomplex, complex_step } ;
@@ -101,71 +105,75 @@ auto get_Ar02(const Model& model, const TType& T, const RhoType& rho, const Mole
}
+template<typename Model, typename Scalar = double, typename VectorType = Eigen::ArrayXd>
+struct VirialDerivatives {
-template <typename Model, typename TType, typename ContainerType>
-typename ContainerType::value_type get_B2vir(const Model& model, const TType T, const ContainerType& molefrac) {
- double h = 1e-100;
- // B_2 = dalphar/drho|T,z at rho=0
- auto B2 = model.alphar(T, std::complex<double>(0.0, h), molefrac).imag()/h;
- return B2;
-}
+ static auto get_B2vir(const Model& model, const Scalar &T, const VectorType& molefrac) {
+ double h = 1e-100;
+ // B_2 = dalphar/drho|T,z at rho=0
+ auto B2 = model.alphar(T, std::complex<double>(0.0, h), molefrac).imag()/h;
+ return B2;
+ }
-/*
-* \f$
-* B_n = \frac{1}{(n-2)!} lim_rho\to 0 d^{n-1}alphar/drho^{n-1}|T,z
-* \f$
-* \param model The model providing the alphar function
-* \param Nderiv The maximum virial coefficient to return; e.g. 5: B_2, B_3, ..., B_5
-* \param T Temperature
-* \param molefrac The mole fractions
-*/
+ /**
+ * \f$
+ * B_n = \frac{1}{(n-2)!} lim_rho\to 0 d^{n-1}alphar/drho^{n-1}|T,z
+ * \f$
+ * \param model The model providing the alphar function
+ * \param Nderiv The maximum virial coefficient to return; e.g. 5: B_2, B_3, ..., B_5
+ * \param T Temperature
+ * \param molefrac The mole fractions
+ */
-template <int Nderiv, ADBackends be = ADBackends::autodiff, typename Model, typename TType, typename ContainerType>
-auto get_Bnvir(const Model& model, const TType &T, const ContainerType& molefrac)
-{
- std::map<int, double> dnalphardrhon;
- if constexpr(be == ADBackends::multicomplex){
- using namespace mcx;
- using fcn_t = std::function<MultiComplex<double>(const MultiComplex<double>&)>;
- fcn_t f = [&model, &T, &molefrac](const auto& rho_) { return model.alphar(T, rho_, molefrac); };
- auto derivs = diff_mcx1(f, 0.0, Nderiv+1, true /* and_val */);
- for (auto n = 1; n <= Nderiv; ++n){
- dnalphardrhon[n] = derivs[n];
+ template <int Nderiv, ADBackends be = ADBackends::autodiff>
+ static auto get_Bnvir(const Model& model, const Scalar &T, const VectorType& molefrac)
+ {
+ std::map<int, double> dnalphardrhon;
+ if constexpr(be == ADBackends::multicomplex){
+ using namespace mcx;
+ using fcn_t = std::function<MultiComplex<double>(const MultiComplex<double>&)>;
+ fcn_t f = [&model, &T, &molefrac](const auto& rho_) { return model.alphar(T, rho_, molefrac); };
+ auto derivs = diff_mcx1(f, 0.0, Nderiv+1, true /* and_val */);
+ for (auto n = 1; n <= Nderiv; ++n){
+ dnalphardrhon[n] = derivs[n];
+ }
}
- }
- else if constexpr(be == ADBackends::autodiff){
- autodiff::HigherOrderDual<Nderiv+1, double> rhodual = 0.0;
- auto f = [&model, &T, &molefrac](const auto& rho_) { return model.alphar(T, rho_, molefrac); };
- auto derivs = derivatives(f, wrt(rhodual), at(rhodual));
- for (auto n = 1; n <= Nderiv; ++n){
- dnalphardrhon[n] = derivs[n];
+ else if constexpr(be == ADBackends::autodiff){
+ autodiff::HigherOrderDual<Nderiv+1, double> rhodual = 0.0;
+ auto f = [&model, &T, &molefrac](const auto& rho_) { return model.alphar(T, rho_, molefrac); };
+ auto derivs = derivatives(f, wrt(rhodual), at(rhodual));
+ for (auto n = 1; n <= Nderiv; ++n){
+ dnalphardrhon[n] = derivs[n];
+ }
}
- }
- else{
- static_assert("algorithmic differentiation backend is invalid");
- }
- std::map<int, TType> o;
- for (int n = 2; n < Nderiv+1; ++n) {
- o[n] = dnalphardrhon[n-1];
- // 0!=1, 1!=1, so only n>3 terms need factorial correction
- if (n > 3) {
- auto factorial = [](int N) {return tgamma(N + 1); };
- o[n] /= factorial(n-2);
+ else{
+ static_assert("algorithmic differentiation backend is invalid");
+ }
+ std::map<int, Scalar> o;
+ for (int n = 2; n < Nderiv+1; ++n) {
+ o[n] = dnalphardrhon[n-1];
+ // 0!=1, 1!=1, so only n>3 terms need factorial correction
+ if (n > 3) {
+ auto factorial = [](int N) {return tgamma(N + 1); };
+ o[n] /= factorial(n-2);
+ }
}
+ return o;
}
- return o;
-}
-template <typename Model, typename TType, typename ContainerType>
-auto get_B12vir(const Model& model, const TType T, const ContainerType& molefrac) {
+ static auto get_B12vir(const Model& model, const Scalar &T, const VectorType& molefrac) {
- auto B2 = get_B2vir(model, T, molefrac); // Overall B2 for mixture
- auto B20 = get_B2vir(model, T, std::valarray<double>({ 1,0 })); // Pure first component with index 0
- auto B21 = get_B2vir(model, T, std::valarray<double>({ 0,1 })); // Pure second component with index 1
- auto z0 = molefrac[0];
- auto B12 = (B2 - z0*z0*B20 - (1-z0)*(1-z0)*B21)/(2*z0*(1-z0));
- return B12;
-}
+ auto B2 = get_B2vir(model, T, molefrac); // Overall B2 for mixture
+ const auto xpure0 = (Eigen::ArrayXd(2) << 1,0).finished();
+ const auto xpure1 = (Eigen::ArrayXd(2) << 0,1).finished();
+ auto B20 = get_B2vir(model, T, xpure0); // Pure first component with index 0
+ auto B21 = get_B2vir(model, T, xpure1); // Pure second component with index 1
+ auto z0 = molefrac[0];
+ auto B12 = (B2 - z0*z0*B20 - (1-z0)*(1-z0)*B21)/(2*z0*(1-z0));
+ return B12;
+ }
+};
+
template<typename Model, typename Scalar = double, typename VectorType = Eigen::ArrayXd>
struct IsochoricDerivatives{
@@ -185,7 +193,7 @@ struct IsochoricDerivatives{
static auto get_pr(const Model& model, const Scalar &T, const VectorType& rhovec)
{
auto rhotot_ = std::accumulate(std::begin(rhovec), std::end(rhovec), (decltype(rhovec[0]))0.0);
- return get_Ar01(model, T, rhotot_, rhovec / rhotot_) * rhotot_ * model.R * T;
+ return ::get_Ar01(model, T, rhotot_, rhovec / rhotot_) * rhotot_ * model.R * T;
}
static auto get_Ar00(const Model& model, const Scalar& T, const VectorType& rhovec) {
diff --git a/interface/pybind11_wrapper.cpp b/interface/pybind11_wrapper.cpp
index ab7bd38..4ba65a1 100644
--- a/interface/pybind11_wrapper.cpp
+++ b/interface/pybind11_wrapper.cpp
@@ -8,9 +8,16 @@
namespace py = pybind11;
+template<typename Model>
+void add_virials(py::module& m) {
+ using vd = VirialDerivatives<Model>;
+ m.def("get_B2vir", &vd::get_B2vir, py::arg("model"), py::arg("T"), py::arg("molefrac"));
+ m.def("get_B12vir", &vd::get_B12vir, py::arg("model"), py::arg("T"), py::arg("molefrac"));
+}
+
template<typename Model>
void add_derivatives(py::module &m) {
- using id = IsochoricDerivatives<Model>;
+ using id = IsochoricDerivatives<Model, double, std::valarray<double>>;
m.def("get_Ar00", &id::get_Ar00, py::arg("model"), py::arg("T"), py::arg("rho"));
m.def("get_Ar10", &id::get_Ar10, py::arg("model"), py::arg("T"), py::arg("rho"));
m.def("get_Psir", &id::get_Psir, py::arg("model"), py::arg("T"), py::arg("rho"));
@@ -20,8 +27,11 @@ void add_derivatives(py::module &m) {
m.def("build_Psir_Hessian_autodiff", &id::build_Psir_Hessian_autodiff, py::arg("model"), py::arg("T"), py::arg("rho"));
m.def("build_Psir_gradient_autodiff", &id::build_Psir_gradient_autodiff, py::arg("model"), py::arg("T"), py::arg("rho"));
+
+ add_virials<Model>(m);
}
+
void init_teqp(py::module& m) {
using vdWEOSd = vdWEOS<double>;
diff --git a/src/multifluid.cpp b/src/multifluid.cpp
index 72a25e1..b180c93 100644
--- a/src/multifluid.cpp
+++ b/src/multifluid.cpp
@@ -83,13 +83,14 @@ int main(){
{
const auto molefrac = (Eigen::ArrayXd(2) << rhovec[0]/rhovec.sum(), rhovec[1]/rhovec.sum()).finished();
- auto B12 = get_B12vir(model, T, molefrac);
+ 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);
const double rho = rhovec.sum();
- volatile double T = 300.0;
+ double T = 300.0;
constexpr int N = 10000;
volatile double alphar;
double rrrr = get_Ar01(model, T, rho, molefrac);
@@ -132,21 +133,21 @@ int main(){
{
Timer t(N);
for (auto i = 0; i < N; ++i) {
- auto o = get_Bnvir<3, ADBackends::autodiff>(model, T, molefrac)[3];
+ auto o = vd::get_Bnvir<3, ADBackends::autodiff>(model, T, molefrac)[3];
}
std::cout << alphar << "; 3 derivs" << std::endl;
}
{
Timer t(N);
for (auto i = 0; i < N; ++i) {
- auto o = get_Bnvir<4, ADBackends::autodiff>(model, T, molefrac)[4];
+ auto o = vd::get_Bnvir<4, ADBackends::autodiff>(model, T, molefrac)[4];
}
std::cout << alphar << "; 4 derivs" << std::endl;
}
{
Timer t(N);
for (auto i = 0; i < N; ++i) {
- auto o = get_Bnvir<5, ADBackends::autodiff>(model, T, molefrac)[5];
+ auto o = vd::get_Bnvir<5, ADBackends::autodiff>(model, T, molefrac)[5];
}
std::cout << alphar << "; 5 derivs" << std::endl;
}
diff --git a/src/tests/catch_tests.cxx b/src/tests/catch_tests.cxx
index 2089c15..fe6a4e9 100644
--- a/src/tests/catch_tests.cxx
+++ b/src/tests/catch_tests.cxx
@@ -44,13 +44,14 @@ TEST_CASE("Check virial coefficients for vdW", "[virial]")
double a = b / ba;
double T = 300.0;
- std::valarray<double> molefrac = { 1.0 };
+ Eigen::ArrayXd molefrac(1); molefrac = 1.0;
constexpr int Nvir = 8;
// Numerical solutions from alphar
- auto Bn = get_Bnvir<Nvir, ADBackends::autodiff>(vdW, T, molefrac);
- auto Bnmcx = get_Bnvir<Nvir, ADBackends::multicomplex>(vdW, T, molefrac);
+ using vd = VirialDerivatives<decltype(vdW)>;
+ auto Bn = vd::get_Bnvir<Nvir, ADBackends::autodiff>(vdW, T, molefrac);
+ auto Bnmcx = vd::get_Bnvir<Nvir, ADBackends::multicomplex>(vdW, T, molefrac);
// Exact solutions for virial coefficients for van der Waals
auto get_vdW_exacts = [a, b, R, T](int Nmax) {
@@ -63,7 +64,7 @@ TEST_CASE("Check virial coefficients for vdW", "[virial]")
auto Bnexact = get_vdW_exacts(Nvir);
// This one with complex step derivatives as another check
- double B2 = get_B2vir(vdW, T, molefrac);
+ double B2 = vd::get_B2vir(vdW, T, molefrac);
double B2exact = b - a / (R * T);
CHECK(std::abs(B2exact-Bnexact[2]) < 1e-15);
CHECK(std::abs(B2-Bnexact[2]) < 1e-15);
@@ -122,7 +123,8 @@ TEST_CASE("Check p three ways for vdW", "[virial][p]")
// Numerical solution from virial expansion
constexpr int Nvir = 8;
- auto Bn = get_Bnvir<Nvir>(model, T, molefrac);
+ using vd = VirialDerivatives<decltype(model)>;
+ auto Bn = vd::get_Bnvir<Nvir>(model, T, molefrac);
auto Z = 1.0;
for (auto i = 2; i <= Nvir; ++i){
Z += Bn[i]*pow(rho, i-1);
@@ -154,8 +156,9 @@ TEST_CASE("Trace critical locus for vdW", "[vdW][crit]")
TEST_CASE("TEST B12", "") {
const auto model = build_vdW();
const double T = 298.15;
- const std::valarray<double> molefrac = { 1/3, 2/3 };
- auto B12 = get_B12vir(model, T, molefrac);
+ const auto molefrac = (Eigen::ArrayXd(2) << 1/3, 2/3).finished();
+ using vd = VirialDerivatives<decltype(model)>;
+ auto B12 = vd::get_B12vir(model, T, molefrac);
}
TEST_CASE("Test psir gradient", "") {
--
GitLab