From 381bd3cc745c621ac5804ef3b4079c4732ba8c28 Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.bell@nist.gov>
Date: Tue, 27 Apr 2021 09:44:22 -0400
Subject: [PATCH] Add script for critical curves in multifluid, VLE routines,
extrapolation from critical
---
CMakeLists.txt | 6 +
include/teqp/algorithms/VLE.hpp | 5 +
include/teqp/algorithms/critical_tracing.hpp | 5 +-
include/teqp/derivs.hpp | 101 ++++---
include/teqp/models/multifluid.hpp | 4 +-
include/teqp/types.hpp | 19 +-
interface/pybind11_wrapper.cpp | 4 +
notebooks/trace_critical.py | 23 +-
src/multifluid.cpp | 301 ++++++++++---------
src/multifluid_crit.cpp | 77 +++++
10 files changed, 344 insertions(+), 201 deletions(-)
create mode 100644 src/multifluid_crit.cpp
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 37f9251..6dcc8dc 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -37,6 +37,12 @@ if (MSVC)
endif()
target_link_libraries(multifluid autodiff)
+add_executable(multifluid2 "${CMAKE_CURRENT_SOURCE_DIR}/src/multifluid2.cpp")
+if (MSVC)
+ target_sources(multifluid2 PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}/externals/Eigen/debug/msvc/eigen.natvis")
+endif()
+target_link_libraries(multifluid2 autodiff)
+
# add_executable(test_autodiff "${CMAKE_CURRENT_SOURCE_DIR}/src/test_autodiff.cpp")
# if (MSVC)
# target_sources(test_autodiff PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}/externals/Eigen/debug/msvc/eigen.natvis")
diff --git a/include/teqp/algorithms/VLE.hpp b/include/teqp/algorithms/VLE.hpp
index c12ddd0..0c94243 100644
--- a/include/teqp/algorithms/VLE.hpp
+++ b/include/teqp/algorithms/VLE.hpp
@@ -101,6 +101,11 @@ auto do_pure_VLE_T(Residual &resid, Scalar rhoL, Scalar rhoV, int maxiter) {
return (Eigen::ArrayXd(2) << rhovec[0], rhovec[1]).finished();
}
+template<typename Model, typename Scalar>
+auto pure_VLE_T(const Model& model, Scalar T, Scalar rhoL, Scalar rhoV, int maxiter) {
+ return do_pure_VLE_T(IsothermPureVLEResiduals(model, T), rhoL, rhoV, maxiter);
+}
+
template<typename Model, typename Scalar>
auto extrapolate_from_critical(const Model& model, const Scalar Tc, const Scalar rhoc, const Scalar T) {
diff --git a/include/teqp/algorithms/critical_tracing.hpp b/include/teqp/algorithms/critical_tracing.hpp
index f59aff2..0186129 100644
--- a/include/teqp/algorithms/critical_tracing.hpp
+++ b/include/teqp/algorithms/critical_tracing.hpp
@@ -51,10 +51,7 @@ struct CriticalTracing {
}
}
- int nonzero_count = 0;
- for (auto& m : mask) {
- nonzero_count += (m == 1);
- }
+ int nonzero_count = mask.count();
auto zero_count = N - nonzero_count;
if (zero_count == 0) {
diff --git a/include/teqp/derivs.hpp b/include/teqp/derivs.hpp
index 1eb6e3f..183ab68 100644
--- a/include/teqp/derivs.hpp
+++ b/include/teqp/derivs.hpp
@@ -12,18 +12,6 @@
#include <autodiff/forward/dual/eigen.hpp>
using namespace autodiff;
-template<typename T>
-auto forceeval(T&& expr)
-{
- using namespace autodiff::detail;
- if constexpr (isDual<T> || isExpr<T> || isNumber<T>) {
- return eval(expr);
- }
- else {
- return expr;
- }
-}
-
/***
* \brief Given a function, use complex step derivatives to calculate the derivative with
* respect to the first variable which here is temperature
@@ -91,7 +79,7 @@ struct TDXDerivatives {
return (1.0/T)*der;
}
else {
- static_assert("algorithmic differentiation backend is invalid in get_Ar10");
+ static_assert(false, "algorithmic differentiation backend is invalid in get_Ar10");
}
}
@@ -116,7 +104,7 @@ struct TDXDerivatives {
return rho*der;
}
else {
- static_assert("algorithmic differentiation backend is invalid in get_Ar01");
+ static_assert(false, "algorithmic differentiation backend is invalid in get_Ar01");
}
}
@@ -129,7 +117,7 @@ struct TDXDerivatives {
return rho*rho*ders[2];
}
else {
- static_assert("algorithmic differentiation backend is invalid in get_Ar02");
+ static_assert(false, "algorithmic differentiation backend is invalid in get_Ar02");
}
}
@@ -141,25 +129,25 @@ struct TDXDerivatives {
auto f = [&model, &T, &molefrac](const auto& rho_) { return eval(model.alphar(T, rho_, molefrac)); };
auto ders = derivatives(f, wrt(rhodual), at(rhodual));
for (auto n = 0; n <= Nderiv; ++n) {
- o[n] = pow(rho, n)*ders[n];
+ o[n] = pow(rho, n) * ders[n];
}
return o;
}
else {
- static_assert("algorithmic differentiation backend is invalid in get_Ar0n");
+ static_assert(false, "algorithmic differentiation backend is invalid in get_Ar0n");
}
}
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)); };
+ 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];
+ return (1 / T) * (1 / T) * ders[2];
}
else {
- static_assert("algorithmic differentiation backend is invalid in get_Ar20");
+ static_assert(false, "algorithmic differentiation backend is invalid in get_Ar20");
}
}
@@ -170,22 +158,22 @@ struct TDXDerivatives {
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);
+ return model.alphar(1.0 / Trecip, rhomolar, molefrac);
};
- std::vector<double> xs = { rho, 1.0/T};
- std::vector<int> order = { 1, 1};
+ 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;
+ return (1.0 / T) * rho * der;
}
else if constexpr (be == ADBackends::autodiff) {
- 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)); };
+ autodiff::dual2nd rhodual = rho, Trecipdual = 1 / T;
+ auto f = [&model, &molefrac](const autodiff::dual2nd& Trecip, const autodiff::dual2nd& 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;
+ return (1.0 / T) * rho * u11;
}
else {
- static_assert("algorithmic differentiation backend is invalid in get_Ar11");
+ static_assert(false, "algorithmic differentiation backend is invalid in get_Ar11");
}
}
@@ -198,20 +186,63 @@ struct TDXDerivatives {
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<double> xs = { rho, 1.0 / T };
std::vector<int> order = { 1, 2 };
auto der = mcx::diff_mcxN(func, xs, order);
- return (1.0/T)*rho*rho*der;
+ return (1.0 / T) * rho * rho * der;
}
else if constexpr (be == ADBackends::autodiff) {
//static_assert("bug in autodiff, can't use autodiff for cross derivative");
autodiff::dual3rd 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 ders = derivatives(f, wrt(Trecipdual, rhodual, rhodual), at(Trecipdual, rhodual));
- return (1.0/T)*rho*rho*ders.back();
+ return (1.0 / T) * rho * rho * ders.back();
+ }
+ else {
+ static_assert(false, "algorithmic differentiation backend is invalid in get_Ar12");
+ }
+ }
+
+ template<ADBackends be = ADBackends::autodiff>
+ static auto get_Ar(const int itau, const int idelta, const Model& model, const Scalar& T, const Scalar& rho, const VectorType& molefrac) {
+ if (itau == 0) {
+ if (idelta == 0) {
+ return get_Ar00(model, T, rho, molefrac);
+ }
+ //else if (idelta == 1) {
+ // return get_Ar01(model, T, rho, molefrac);
+ //}
+ //else if (idelta == 2) {
+ // return get_Ar02(model, T, rho, molefrac);
+ //}
+ else {
+ throw std::invalid_argument("Invalid value for idelta");
+ }
+ }
+ else if (itau == 1){
+ if (idelta == 0) {
+ return get_Ar10(model, T, rho, molefrac);
+ }
+ //else if (idelta == 1) {
+ // return get_Ar11(model, T, rho, molefrac);
+ //}
+ /*else if (idelta == 2) {
+ return get_Ar12(model, T, rho, molefrac);
+ }*/
+ else {
+ throw std::invalid_argument("Invalid value for idelta");
+ }
+ }
+ else if (itau == 2) {
+ if (idelta == 0) {
+ return get_Ar20(model, T, rho, molefrac);
+ }
+ else {
+ throw std::invalid_argument("Invalid value for idelta");
+ }
}
else {
- static_assert("algorithmic differentiation backend is invalid in get_Ar12");
+ throw std::invalid_argument("Invalid value for idelta");
}
}
@@ -219,7 +250,7 @@ struct TDXDerivatives {
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);
+ auto Ar20 = get_Ar20<be>(model, T, rho, molefrac);
return -3*(Ar01-Ar11)/Ar20;
}
};
@@ -266,7 +297,7 @@ struct VirialDerivatives {
}
}
else{
- static_assert("algorithmic differentiation backend is invalid");
+ static_assert(false, "algorithmic differentiation backend is invalid");
}
std::map<int, Scalar> o;
for (int n = 2; n < Nderiv+1; ++n) {
diff --git a/include/teqp/models/multifluid.hpp b/include/teqp/models/multifluid.hpp
index 5c1626a..328d90f 100644
--- a/include/teqp/models/multifluid.hpp
+++ b/include/teqp/models/multifluid.hpp
@@ -6,6 +6,8 @@
#include <fstream>
#include <string>
#include <cmath>
+#include <optional>
+#include "teqp/types.hpp"
#include "MultiComplex/MultiComplex.hpp"
// See https://eigen.tuxfamily.org/dox/TopicCustomizing_CustomScalar.html
@@ -75,7 +77,7 @@ public:
const CorrespondingTerm corr;
const DepartureTerm dep;
- const double R = get_R_gas<double>();
+ const double R = 1.380649e-23 * 6.02214076e23; ///< Exact value, given by k_B*N_A
MultiFluid(ReducingFunction&& redfunc, CorrespondingTerm&& corr, DepartureTerm&& dep) : redfunc(redfunc), corr(corr), dep(dep) {};
diff --git a/include/teqp/types.hpp b/include/teqp/types.hpp
index 91a54a0..9061d31 100644
--- a/include/teqp/types.hpp
+++ b/include/teqp/types.hpp
@@ -9,4 +9,21 @@ template <typename Container>
struct is_container : std::false_type { };
template <typename... Ts> struct is_container<std::vector<Ts...> > : std::true_type { };
template <typename... Ts> struct is_container<std::valarray<Ts...> > : std::true_type { };
-// Missing Eigen::Array and Eigen::Matrix types here
\ No newline at end of file
+// Missing Eigen::Array and Eigen::Matrix types here
+
+// autodiff include
+#include <autodiff/forward/dual.hpp>
+#include <autodiff/forward/dual/eigen.hpp>
+using namespace autodiff;
+
+template<typename T>
+auto forceeval(T&& expr)
+{
+ using namespace autodiff::detail;
+ if constexpr (isDual<T> || isExpr<T> || isNumber<T>) {
+ return eval(expr);
+ }
+ else {
+ return expr;
+ }
+}
\ No newline at end of file
diff --git a/interface/pybind11_wrapper.cpp b/interface/pybind11_wrapper.cpp
index 1680037..bb368bc 100644
--- a/interface/pybind11_wrapper.cpp
+++ b/interface/pybind11_wrapper.cpp
@@ -13,6 +13,7 @@
#include "teqp/algorithms/critical_tracing.hpp"
#include "teqp/models/pcsaft.hpp"
#include "teqp/models/multifluid.hpp"
+#include "teqp/algorithms/VLE.hpp"
namespace py = pybind11;
@@ -50,6 +51,9 @@ void add_derivatives(py::module &m, Wrapper &cls) {
using ct = CriticalTracing<Model, double, Eigen::Array<double, Eigen::Dynamic, 1>>;
m.def("trace_critical_arclength_binary", &ct::trace_critical_arclength_binary);
+ m.def("extrapolate_from_critical", &extrapolate_from_critical<Model, double>);
+ m.def("pure_VLE_T", &pure_VLE_T<Model, double>);
+
//cls.def("get_Ar01", [](const Model& m, const double T, const Eigen::ArrayXd& rhovec) { return id::get_Ar01(m, T, rhovec); });
//cls.def("get_Ar10", [](const Model& m, const double T, const Eigen::ArrayXd& rhovec) { return id::get_Ar10(m, T, rhovec); });
using tdx = TDXDerivatives<Model, double, Eigen::Array<double, Eigen::Dynamic, 1> >;
diff --git a/notebooks/trace_critical.py b/notebooks/trace_critical.py
index 10e1881..a953e18 100644
--- a/notebooks/trace_critical.py
+++ b/notebooks/trace_critical.py
@@ -9,7 +9,7 @@ import pandas
import scipy.optimize
def trace():
- model = teqp.build_multifluid_model(["Methane", "Ethane"], '../mycp', '../mycp/dev/mixtures/mixture_binary_pairs.json')
+ model = teqp.build_multifluid_model(["R1234yf", "R32"], '../mycp', '../mycp/dev/mixtures/mixture_binary_pairs.json')
Tcvec = model.get_Tcvec()
rhocvec = 1/model.get_vcvec()
models = [model]
@@ -18,7 +18,7 @@ def trace():
for i in range(2):
z = [0,0]; z[i] = 1
z = np.array(z)
- p = 8.314462618*Tcvec[i]*rhocvec[i]*(1+model.get_Ar01(Tcvec[i], rhocvec[i]*z))
+ p = 8.314462618*Tcvec[i]*rhocvec[i]*(1+model.get_Ar01(Tcvec[i], rhocvec[i], z))
pc_Pa.append(p)
model = teqp.vdWEOS(Tcvec, pc_Pa)
models.append(model)
@@ -26,24 +26,27 @@ def trace():
i = 1
T0 = Tcvec[i]
z = np.array([0, 0]); z[i] = 1
- if isinstance(model, teqp.vdWEOS):
- rho0 = z*pc_Pa[i]/(8.314462618*Tcvec[i])*8/3
- else:
- rho0 = z*rhocvec[i]
fig1, ax1 = plt.subplots(1,1)
fig2, ax2 = plt.subplots(1,1)
fig3, ax3 = plt.subplots(1,1)
for model in models:
+
+ if isinstance(model, teqp.vdWEOS):
+ rho0 = z*pc_Pa[i]/(8.314462618*Tcvec[i])*8/3
+ else:
+ rho0 = z*rhocvec[i]
+
tic = timeit.default_timer()
curveJSON = teqp.trace_critical_arclength_binary(model, T0, rho0, "")
toc = timeit.default_timer()
- print(toc-tic, 's')
+ print(toc-tic)
+
df = pandas.DataFrame(curveJSON)
df['z0 / mole frac.'] = df['rho0 / mol/m^3']/(df['rho0 / mol/m^3']+df['rho1 / mol/m^3'])
- ax1.plot(df['z0 / mole frac.'], df['T / K'])
- ax2.plot(df['T / K'], df['p / Pa'])
- ax3.plot(df['z0 / mole frac.'], df['rho0 / mol/m^3']+df['rho1 / mol/m^3'])
+ ax1.plot(df['z0 / mole frac.'], df['T / K']); ax1.set(xlabel='z0 / mole frac.', ylabel='$T$ / K')
+ ax2.plot(df['T / K'], df['p / Pa']); ax2.set(xlabel='$T$ / K', ylabel='$p$ / Pa')
+ ax3.plot(df['z0 / mole frac.'], df['s^+']/df['s^+'].iloc[0]); ax3.set(xlabel='x0 / mole frac.', ylabel='')
plt.show()
if __name__ == '__main__':
diff --git a/src/multifluid.cpp b/src/multifluid.cpp
index 44461a1..cbf7919 100644
--- a/src/multifluid.cpp
+++ b/src/multifluid.cpp
@@ -5,146 +5,146 @@
#include "teqp/algorithms/critical_tracing.hpp"
#include <optional>
-
-class Timer {
-private:
- int N;
- decltype(std::chrono::steady_clock::now()) tic;
-public:
- Timer(int N) : N(N), tic(std::chrono::steady_clock::now()){}
- ~Timer() {
- auto elap = std::chrono::duration<double>(std::chrono::steady_clock::now()-tic).count();
- std::cout << elap/N*1e6 << " us/call" << std::endl;
- }
-};
-
-void trace_critical_loci(const std::string &coolprop_root, const nlohmann::json &BIPcollection) {
- std::vector<std::vector<std::string>> pairs = {
- { "CarbonDioxide", "R1234YF" }, { "CarbonDioxide","R1234ZE(E)" }, { "ETHYLENE","R1243ZF" },
- { "R1234YF","R1234ZE(E)" }, { "R134A","R1234YF" }, { "R23","R1234YF" },
- { "R32","R1123" }, { "R32","R1234YF" }, { "R32","R1234ZE(E)" }
- };
- for (auto &pp : pairs) {
- using ModelType = decltype(build_multifluid_model(pp, coolprop_root, BIPcollection));
- std::optional<ModelType> optmodel{std::nullopt};
- try {
- optmodel.emplace(build_multifluid_model(pp, coolprop_root, BIPcollection));
- }
- catch (std::exception &e) {
- std::cout << e.what() << std::endl;
- std::cout << pp[0] << "&" << pp[1] << std::endl;
- continue;
- }
- for (int i : {0, 1}){
- const auto &model = optmodel.value();
- auto rhoc0 = 1.0 / model.redfunc.vc[i];
- auto T0 = model.redfunc.Tc[i];
- Eigen::ArrayXd rhovec(2); rhovec[i] = { rhoc0 }; rhovec[1L - i] = 0.0;
-
- using ct = CriticalTracing<ModelType>;
-
- // Non-analytic terms make it impossible to initialize AT the pure components
- if (pp[0] == "CarbonDioxide" || pp[1] == "CarbonDioxide") {
- if (i == 0) {
- rhovec[i] *= 0.9999;
- rhovec[1L - i] = 0.9999;
- }
- else {
- rhovec[i] *= 1.0001;
- rhovec[1L - i] = 1.0001;
- }
- double zi = rhovec[i] / rhovec.sum();
- double T = zi * model.redfunc.Tc[i] + (1 - zi) * model.redfunc.Tc[1L - i];
- double z0 = (i == 0) ? zi : 1-zi;
- auto [Tnew, rhonew] = ct::critical_polish_molefrac(model, T, rhovec, z0);
- T0 = Tnew;
- rhoc0 = rhovec.sum();
- }
- std::string filename = pp[0] + "_" + pp[1] + ".csv";
- ct::trace_critical_arclength_binary(model, T0, rhovec, filename);
- }
- }
-}
-
-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;
- double T = 300;
- {
- 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);
-
- const double rho = rhovec.sum();
- double T = 300.0;
- constexpr int N = 10000;
- volatile double alphar;
- 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) {
- alphar = model.alphar(T, rho, molefrac);
- }
- std::cout << alphar << " function call" << std::endl;
- }
- {
- Timer t(N);
- for (auto i = 0; i < N; ++i) {
- 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 = 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 = 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 = tdx::get_Ar02(model, T, rho, molefrac);
- }
- std::cout << alphar << "; 2nd autodiff" << std::endl;
- }
- {
- Timer t(N);
- for (auto i = 0; i < N; ++i) {
- auto o = vd::template 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 = vd::template 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 = vd::template get_Bnvir<5, ADBackends::autodiff>(model, T, molefrac)[5];
- }
- std::cout << alphar << "; 5 derivs" << std::endl;
- }
- }
-}
+//
+//class Timer {
+//private:
+// int N;
+// decltype(std::chrono::steady_clock::now()) tic;
+//public:
+// Timer(int N) : N(N), tic(std::chrono::steady_clock::now()){}
+// ~Timer() {
+// auto elap = std::chrono::duration<double>(std::chrono::steady_clock::now()-tic).count();
+// std::cout << elap/N*1e6 << " us/call" << std::endl;
+// }
+//};
+//
+//void trace_critical_loci(const std::string &coolprop_root, const nlohmann::json &BIPcollection) {
+// std::vector<std::vector<std::string>> pairs = {
+// { "CarbonDioxide", "R1234YF" }, { "CarbonDioxide","R1234ZE(E)" }, { "ETHYLENE","R1243ZF" },
+// { "R1234YF","R1234ZE(E)" }, { "R134A","R1234YF" }, { "R23","R1234YF" },
+// { "R32","R1123" }, { "R32","R1234YF" }, { "R32","R1234ZE(E)" }
+// };
+// for (auto &pp : pairs) {
+// using ModelType = decltype(build_multifluid_model(pp, coolprop_root, BIPcollection));
+// std::optional<ModelType> optmodel{std::nullopt};
+// try {
+// optmodel.emplace(build_multifluid_model(pp, coolprop_root, BIPcollection));
+// }
+// catch (std::exception &e) {
+// std::cout << e.what() << std::endl;
+// std::cout << pp[0] << "&" << pp[1] << std::endl;
+// continue;
+// }
+// for (int i : {0, 1}){
+// const auto &model = optmodel.value();
+// auto rhoc0 = 1.0 / model.redfunc.vc[i];
+// auto T0 = model.redfunc.Tc[i];
+// Eigen::ArrayXd rhovec(2); rhovec[i] = { rhoc0 }; rhovec[1L - i] = 0.0;
+//
+// using ct = CriticalTracing<ModelType>;
+//
+// // Non-analytic terms make it impossible to initialize AT the pure components
+// if (pp[0] == "CarbonDioxide" || pp[1] == "CarbonDioxide") {
+// if (i == 0) {
+// rhovec[i] *= 0.9999;
+// rhovec[1L - i] = 0.9999;
+// }
+// else {
+// rhovec[i] *= 1.0001;
+// rhovec[1L - i] = 1.0001;
+// }
+// double zi = rhovec[i] / rhovec.sum();
+// double T = zi * model.redfunc.Tc[i] + (1 - zi) * model.redfunc.Tc[1L - i];
+// double z0 = (i == 0) ? zi : 1-zi;
+// auto [Tnew, rhonew] = ct::critical_polish_molefrac(model, T, rhovec, z0);
+// T0 = Tnew;
+// rhoc0 = rhovec.sum();
+// }
+// std::string filename = pp[0] + "_" + pp[1] + ".csv";
+// ct::trace_critical_arclength_binary(model, T0, rhovec, filename);
+// }
+// }
+//}
+//
+//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;
+// double T = 300;
+// {
+// 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);
+//
+// const double rho = rhovec.sum();
+// double T = 300.0;
+// constexpr int N = 10000;
+// volatile double alphar;
+// 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) {
+// alphar = model.alphar(T, rho, molefrac);
+// }
+// std::cout << alphar << " function call" << std::endl;
+// }
+// {
+// Timer t(N);
+// for (auto i = 0; i < N; ++i) {
+// 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 = 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 = 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 = tdx::get_Ar02(model, T, rho, molefrac);
+// }
+// std::cout << alphar << "; 2nd autodiff" << std::endl;
+// }
+// {
+// Timer t(N);
+// for (auto i = 0; i < N; ++i) {
+// auto o = vd::template 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 = vd::template 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 = vd::template get_Bnvir<5, ADBackends::autodiff>(model, T, molefrac)[5];
+// }
+// std::cout << alphar << "; 5 derivs" << std::endl;
+// }
+// }
+//}
int main(){
@@ -153,10 +153,10 @@ int main(){
auto BIPcollection = coolprop_root + "/dev/mixtures/mixture_binary_pairs.json";
// Critical curves
- {
+ /* {
Timer t(1);
trace_critical_loci(coolprop_root, BIPcollection);
- }
+ }*/
//time_calls(coolprop_root, BIPcollection);
@@ -166,18 +166,19 @@ int main(){
double T = 300;
const auto molefrac = rhovec/rhovec.sum();
- using tdx = TDXDerivatives<decltype(model)>;
+ using tdx = TDXDerivatives<decltype(model), double, Eigen::ArrayXd>;
const auto b = ADBackends::autodiff;
- auto rho = rhovec.sum();
+ auto rho = rhovec.sum();
auto alphar = model.alphar(T, rho, rhovec);
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 Ar10 = tdx::get_Ar10<b>(model, T, rho, molefrac);
+ auto Ar02 = tdx::get_Ar02<b>(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);
+
+ //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), double, Eigen::ArrayXd>;
+ //auto splus = id::get_splus(model, T, rhovec);*/
int ttt = 0;
}
diff --git a/src/multifluid_crit.cpp b/src/multifluid_crit.cpp
new file mode 100644
index 0000000..4823cd3
--- /dev/null
+++ b/src/multifluid_crit.cpp
@@ -0,0 +1,77 @@
+#define USE_AUTODIFF
+
+#include "teqp/core.hpp"
+#include "teqp/models/multifluid.hpp"
+#include "teqp/algorithms/critical_tracing.hpp"
+
+#include <optional>
+
+class Timer {
+private:
+ int N;
+ decltype(std::chrono::steady_clock::now()) tic;
+public:
+ Timer(int N) : N(N), tic(std::chrono::steady_clock::now()){}
+ ~Timer() {
+ auto elap = std::chrono::duration<double>(std::chrono::steady_clock::now()-tic).count();
+ std::cout << elap/N*1e6 << " us/call" << std::endl;
+ }
+};
+
+void trace_critical_loci(const std::string &coolprop_root, const nlohmann::json &BIPcollection) {
+ std::vector<std::vector<std::string>> pairs = {
+ { "CarbonDioxide", "R1234YF" }, { "CarbonDioxide","R1234ZE(E)" }, { "ETHYLENE","R1243ZF" },
+ { "R1234YF","R1234ZE(E)" }, { "R134A","R1234YF" }, { "R23","R1234YF" },
+ { "R32","R1123" }, { "R32","R1234YF" }, { "R32","R1234ZE(E)" }
+ };
+ for (auto &pp : pairs) {
+ using ModelType = decltype(build_multifluid_model(pp, coolprop_root, BIPcollection));
+ std::optional<ModelType> optmodel{std::nullopt};
+ try {
+ optmodel.emplace(build_multifluid_model(pp, coolprop_root, BIPcollection));
+ }
+ catch (std::exception &e) {
+ std::cout << e.what() << std::endl;
+ std::cout << pp[0] << "&" << pp[1] << std::endl;
+ continue;
+ }
+ for (int i : {0, 1}){
+ const auto &model = optmodel.value();
+ auto rhoc0 = 1.0 / model.redfunc.vc[i];
+ auto T0 = model.redfunc.Tc[i];
+ Eigen::ArrayXd rhovec(2); rhovec[i] = { rhoc0 }; rhovec[1L - i] = 0.0;
+
+ using ct = CriticalTracing<ModelType>;
+
+ // Non-analytic terms make it impossible to initialize AT the pure components
+ if (pp[0] == "CarbonDioxide" || pp[1] == "CarbonDioxide") {
+ if (i == 0) {
+ rhovec[i] *= 0.9999;
+ rhovec[1L - i] = 0.9999;
+ }
+ else {
+ rhovec[i] *= 1.0001;
+ rhovec[1L - i] = 1.0001;
+ }
+ double zi = rhovec[i] / rhovec.sum();
+ double T = zi * model.redfunc.Tc[i] + (1 - zi) * model.redfunc.Tc[1L - i];
+ double z0 = (i == 0) ? zi : 1-zi;
+ auto [Tnew, rhonew] = ct::critical_polish_molefrac(model, T, rhovec, z0);
+ T0 = Tnew;
+ rhoc0 = rhovec.sum();
+ }
+ std::string filename = pp[0] + "_" + pp[1] + ".csv";
+ ct::trace_critical_arclength_binary(model, T0, rhovec, filename);
+ }
+ }
+}
+
+int main(){
+
+ std::string coolprop_root = "../mycp";
+ auto BIPcollection = coolprop_root + "/dev/mixtures/mixture_binary_pairs.json";
+ // Critical curves
+ Timer t(1);
+ trace_critical_loci(coolprop_root, BIPcollection);
+ return EXIT_SUCCESS;
+}
\ No newline at end of file
--
GitLab