From f0c132b9433a02476b3cfa2fa83627c980029866 Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.bell@nist.gov>
Date: Wed, 24 Mar 2021 16:07:04 -0400
Subject: [PATCH] Correct the virial coefficient calculation
---
include/teqp/derivs.hpp | 19 +++++++++++++++++--
src/main.cpp | 18 +++++++++---------
2 files changed, 26 insertions(+), 11 deletions(-)
diff --git a/include/teqp/derivs.hpp b/include/teqp/derivs.hpp
index 51dd950..6d28070 100644
--- a/include/teqp/derivs.hpp
+++ b/include/teqp/derivs.hpp
@@ -3,6 +3,7 @@
#include <optional>
#include <complex>
#include <tuple>
+#include <map>
#include "MultiComplex/MultiComplex.hpp"
@@ -110,13 +111,27 @@ typename ContainerType::value_type get_B2vir(const Model& model, const TType T,
return B2;
}
+/*
+ B_n = \frac{1}{(n-2)!} lim_rho\to 0 d^{n-1}alphar/drho^{n-1}|T,z
+*/
+
template <typename Model, typename TType, typename ContainerType>
auto get_Bnvir(const Model& model, int Nderiv, const TType T, const ContainerType& molefrac) {
- // B_n = lim_rho\to 0 d^{n-1}alphar/drho^{n-1}|T,z
+
using namespace mcx;
using fcn_t = std::function<MultiComplex<double>(const MultiComplex<double>&)>;
fcn_t f = [&model, &T, &molefrac](const MultiComplex<double>& rho_) -> MultiComplex<double> { return model.alphar(T, rho_, molefrac); };
- return diff_mcx1(f, 0.0, Nderiv, true /* and_val */);
+ std::map<int, TType> o;
+ auto dalphardrhon = diff_mcx1(f, 0.0, Nderiv+1, true /* and_val */);
+ for (int n = 2; n < Nderiv+1; ++n) {
+ o[n] = dalphardrhon[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;
}
/***
diff --git a/src/main.cpp b/src/main.cpp
index e7447cb..1ef6bba 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -42,17 +42,17 @@ void test_vdW() {
auto Nvir = 8;
auto Bn = get_Bnvir(vdW, Nvir, T, molefrac);
// Exact solutions for virial coefficients for van der Waals
- auto get_vdW_exacts = [a,b,R,T](int additional){
- std::vector<double> o = {0, b - a / (R * T)};
- for (auto i = 0; i < additional; ++i) {
- auto factorial = [](int N){return tgamma(N + 1);};
- o.push_back(factorial(i+1)*pow(b,i+2));
+ auto get_vdW_exacts = [a,b,R,T](int Nmax){
+ std::map<int, double> o = {{2, b - a / (R * T)}};
+ for (auto i = 3; i <= Nmax; ++i) {
+ o[i] = pow(b, i-1);
}
- return std::valarray<double>(&(o[0]), o.size());
+ return o;
};
- std::valarray<double> Bnexact = get_vdW_exacts(Nvir-1);
- std::valarray<double> errrr = Bnexact - std::valarray<double>(&(Bn[0]), Bn.size());
-
+ auto Bnexact = get_vdW_exacts(Nvir);
+ for (auto i = 2; i <= Nvir; ++i){
+ std::cout << std::scientific << i << ", " << Bnexact[i] << ", " << Bn[i] << ", " << std::abs(Bnexact[i]-Bn[i]) << std::endl;
+ }
int rr = 0;
}
--
GitLab