diff --git a/include/teqp/models/multifluid.hpp b/include/teqp/models/multifluid.hpp
index 7b74a0bb63a5bf8e653d597b21f054537a25fbc8..8e8685b7dfd10b941d5f34af3a66f943096ce3a5 100644
--- a/include/teqp/models/multifluid.hpp
+++ b/include/teqp/models/multifluid.hpp
@@ -51,6 +51,16 @@ public:
}
return alphar;
}
+
+ template<typename TauType, typename DeltaType>
+ auto alphari(const TauType& tau, const DeltaType& delta, std::size_t i) const {
+ using resulttype = std::common_type_t<decltype(tau), decltype(delta)>; // Type promotion, without the const-ness
+ return EOSs[i].alphar(tau, delta);
+ }
+
+ auto get_EOS(std::size_t i) const{
+ return EOSs[i];
+ }
};
template<typename FCollection, typename DepartureFunctionCollection>
diff --git a/include/teqp/models/multifluid_eosterms.hpp b/include/teqp/models/multifluid_eosterms.hpp
index c9cf115616c9a94eb074eaaa27f30763498702c1..83714d81d4acebbd75c694af506ae478539421cc 100644
--- a/include/teqp/models/multifluid_eosterms.hpp
+++ b/include/teqp/models/multifluid_eosterms.hpp
@@ -7,7 +7,12 @@ public:
template<typename TauType, typename DeltaType>
auto alphar(const TauType& tau, const DeltaType& delta) const {
- return forceeval((n * exp(t * log(tau) + d * log(delta) - c * powIVi(delta, l_i))).sum());
+ using result = std::common_type_t<TauType, DeltaType>;
+ result r = 0.0, lntau = log(tau), lndelta = log(delta);
+ for (auto i = 0; i < n.size(); ++i) {
+ r += n[i] * exp(t[i]*lntau + d[i]*lndelta -c[i]*powi(delta, l_i[i]));
+ }
+ return r;
}
};
diff --git a/include/teqp/types.hpp b/include/teqp/types.hpp
index 0c85275bb815596cbcb412cfee82c1ce8b6222f4..47505dfb53bb12c8eb2d4f0ff3d2528eda7c92c2 100644
--- a/include/teqp/types.hpp
+++ b/include/teqp/types.hpp
@@ -83,9 +83,20 @@ public:
/// From Ulrich Deiters
template <typename T> // arbitrary integer power
T powi(const T& x, int n) {
- if (n == 0)
+ switch (n) {
+ case 0:
return static_cast<T>(1.0); // x^0 = 1 even for x == 0
- else if (n < 0){
+ case 1:
+ return static_cast<T>(x);
+ case 2:
+ return static_cast<T>(x*x);
+ case 3:
+ return static_cast<T>(x*x*x);
+ case 4:
+ auto x2 = x * x;
+ return static_cast<T>(x2*x2);
+ }
+ if (n < 0){
using namespace autodiff::detail;
if constexpr (isDual<T> || isExpr<T>) {
return eval(powi(eval(1.0/x), -n));
diff --git a/src/time_REFPROP.cpp b/src/time_REFPROP.cpp
index a499a527c8cbfb8d44602815005706f3e1bf496e..01e73a562c93e5d9e93a53bf1e2dbc4f38d39a24 100644
--- a/src/time_REFPROP.cpp
+++ b/src/time_REFPROP.cpp
@@ -13,6 +13,8 @@
#include <random>
#include <numeric>
+#include "teqp/models/multifluid.hpp"
+
int main()
{
// You may need to change this path to suit your installation
@@ -27,8 +29,14 @@ int main()
if (!loaded_REFPROP){return EXIT_FAILURE; }
SETPATHdll(const_cast<char*>(path.c_str()), 400);
int ierr = 0, nc = 1;
- char herr[255], hfld[10000] = "METHANE", hhmx[255] = "HMX.BNC", href[4] = "DEF";
+ char herr[255], hfld[10000] = "KRYPTON", hhmx[255] = "HMX.BNC", href[4] = "DEF";
SETUPdll(nc, hfld, hhmx, href, ierr, herr, 10000, 255, 3, 255);
+ {
+ char hflag[256] = "Cache ";
+ int jFlag = 3, kFlag = -1;
+ FLAGSdll(hflag, jFlag, kFlag, ierr, herr, 255, 255);
+ std::cout << kFlag << std::endl;
+ }
if (ierr != 0) printf("This ierr: %d herr: %s\n", ierr, herr);
{
@@ -37,25 +45,52 @@ int main()
// Random speed testing
std::default_random_engine re;
- std::valarray<double> Ts(100000);
+ std::valarray<double> taus(100000);
{
- std::uniform_real_distribution<double> unif(89, 360);
- std::transform(std::begin(Ts), std::end(Ts), std::begin(Ts), [&unif, &re](double x) { return unif(re); });
+ std::uniform_real_distribution<double> unif(2.0941098901098902, 2.1941098901098902);
+ std::transform(std::begin(taus), std::end(taus), std::begin(taus), [&unif, &re](double x) { return unif(re); });
}
- std::valarray<double> ds(Ts.size()); {
- std::uniform_real_distribution<double> unif(1, 1.001);
- std::transform(std::begin(ds), std::end(ds), std::begin(ds), [&unif, &re](double x) { return unif(re); });
+ std::valarray<double> deltas(taus.size()); {
+ std::uniform_real_distribution<double> unif(0.0015981745536338204, 0.0016981745536338204);
+ std::transform(std::begin(deltas), std::end(deltas), std::begin(deltas), [&unif, &re](double x) { return unif(re); });
}
int itau = 0, idelta = 0;
- for (auto repeat = 0; repeat < 10; ++repeat) {
- std::valarray<double> ps = 0.0 * Ts;
+ for (auto repeat = 0; repeat < 100; ++repeat) {
+ std::valarray<double> ps = 0.0 * taus;
+ double rhoc = 10139.128;
+ double Tc = 190.564;
auto tic = std::chrono::high_resolution_clock::now();
- for (auto i = 0; i < Ts.size(); ++i) {
- PHIXdll(itau, idelta, Ts[i], ds[i], z, p); ps[i] = p;
+ for (auto i = 0; i < taus.size(); ++i) {
+ double Ar01 = -10000;
+ PHIXdll(itau, idelta, taus[i], deltas[i], z, Ar01); ps[i] = Ar01; //ps[i] = (Ar01 + 1)*deltas[i]*rhoc*8.31451*Tc/taus[i];
+ /*idelta = 1;
+ PHIXdll(itau, idelta, Ts[i], ds[i], z, p); ps[i] += p;
+ idelta = 2;
+ PHIXdll(itau, idelta, Ts[i], ds[i], z, p); ps[i] += p;
+ idelta = 3;
+ PHIXdll(itau, idelta, Ts[i], ds[i], z, p); ps[i] += p;*/
+ //PRESSdll(Ts[i], ds[i], z, p); ps[i] = p;
}
auto toc = std::chrono::high_resolution_clock::now();
- double elap_us = std::chrono::duration<double>(toc - tic).count() / Ts.size() * 1e6;
- std::cout << elap_us << " us/call for PRESSdll; " << std::accumulate(std::begin(ps), std::end(ps), 0.0) << std::endl;
+ double elap_us = std::chrono::duration<double>(toc - tic).count() / taus.size() * 1e6;
+ std::cout << "Repeat " << std::to_string(repeat) << ": " << elap_us << " us/call for PRESSdll; " << std::accumulate(std::begin(ps), std::end(ps), 0.0)/ps.size() << std::endl;
+ }
+
+ // And the same example with teqp
+ auto model = build_multifluid_model({ "Krypton" }, "../mycp", "../mycp/dev/mixtures/mixture_binary_pairs.json");
+ auto N = taus.size();
+ auto c = (Eigen::ArrayXd(1) << 1.0).finished();
+
+ // Get reference to pure fluid model
+ const auto f = model.corr.get_EOS(0);
+ for (auto counter = 0; counter < 100; ++counter)
+ {
+ double o = 0.0;
+ Timer t(N);
+ for (auto j = 0; j < N; ++j) {
+ o += f.alphar(taus[j], deltas[j]);
+ }
+ std::cout << o/N << std::endl;
}
}
return EXIT_SUCCESS;