From 9a518533d14b1835bd015b45e7b4a51817257028 Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.bell@nist.gov>
Date: Fri, 18 Jun 2021 15:11:05 -0400
Subject: [PATCH] Departure functions refactored and fixed
---
include/teqp/models/multifluid.hpp | 235 ++++++++++----------
include/teqp/models/multifluid_eosterms.hpp | 59 +++--
include/teqp/types.hpp | 22 +-
3 files changed, 177 insertions(+), 139 deletions(-)
diff --git a/include/teqp/models/multifluid.hpp b/include/teqp/models/multifluid.hpp
index 40441dc..953fc11 100644
--- a/include/teqp/models/multifluid.hpp
+++ b/include/teqp/models/multifluid.hpp
@@ -8,7 +8,6 @@
#include <cmath>
#include <optional>
#include <variant>
-#include <set>
#include "teqp/types.hpp"
#include "MultiComplex/MultiComplex.hpp"
@@ -268,114 +267,137 @@ public:
template<typename MoleFractions> auto get_rhor(const MoleFractions& molefracs) const { return 1.0 / Y(molefracs, vc, betaV, Yv); }
};
-class MultiFluidDepartureFunction {
-public:
- enum class types { NOTSETTYPE, GERG2004, GaussianExponential, NoDeparture };
-private:
- types type = types::NOTSETTYPE;
-public:
- Eigen::ArrayXd n, t, d, c, l, eta, beta, gamma, epsilon;
+inline auto get_departure_function_matrix(const std::string& coolprop_root, const nlohmann::json& BIPcollection, const std::vector<std::string>& components, const nlohmann::json& flags) {
- void set_type(const std::string& kind) {
- if (kind == "GERG-2004" || kind == "GERG-2008") {
- type = types::GERG2004;
- }
- else if (kind == "Gaussian+Exponential") {
- type = types::GaussianExponential;
+ // Allocate the matrix with default models
+ std::vector<std::vector<DepartureTerms>> funcs(components.size()); for (auto i = 0; i < funcs.size(); ++i) { funcs[i].resize(funcs.size()); }
+
+ // Load the collection of data on departure functions
+ auto depcollection = nlohmann::json::parse(std::ifstream(coolprop_root + "/dev/mixtures/mixture_departure_functions.json"));
+ auto get_departure_json = [&depcollection](const std::string& Name) {
+ for (auto& el : depcollection) {
+ if (el["Name"] == Name) { return el; }
}
- else if (kind == "none") {
- type = types::NoDeparture;
+ throw std::invalid_argument("Bad argument");
+ };
+
+ auto build_power = [&](auto term) {
+ std::size_t N = term["n"].size();
+
+ PowerEOSTerm eos;
+
+ auto eigorzero = [&term, &N](const std::string& name) -> Eigen::ArrayXd {
+ if (!term[name].empty()) {
+ return toeig(term[name]);
+ }
+ else {
+ return Eigen::ArrayXd::Zero(N);
+ }
+ };
+
+
+ eos.n = eigorzero("n");
+ eos.t = eigorzero("t");
+ eos.d = eigorzero("d");
+
+ Eigen::ArrayXd c(N), l(N); c.setZero();
+ if (term["l"].empty()) {
+ // exponential part not included
+ l.setZero();
}
else {
- throw std::invalid_argument("Bad type:" + kind);
- }
- }
-
- template<typename TauType, typename DeltaType>
- auto alphar(const TauType& tau, const DeltaType& delta) const {
- switch (type) {
- case (types::GaussianExponential):
- return forceeval((n * exp(t*log(tau) + d*log(delta)-c*pow(delta, l)-eta * (delta - epsilon).square() - beta * (tau - gamma).square())).sum());
- case (types::GERG2004):
- return forceeval((n * exp(t*log(tau) + d*log(delta) -eta * (delta - epsilon).square() - beta * (delta - gamma))).sum());
- case (types::NoDeparture):
- return forceeval(0.0*(tau*delta));
- default:
- throw - 1;
+ l = toeig(term["l"]);
+ // l is included, use it to build c; c_i = 1 if l_i > 0, zero otherwise
+ for (auto i = 0; i < c.size(); ++i) {
+ if (l[i] > 0) {
+ c[i] = 1.0;
+ }
+ }
}
- }
-};
+ eos.c = c;
+ eos.l = l;
-inline auto get_departure_function_matrix(const std::string& coolprop_root, const nlohmann::json& BIPcollection, const std::vector<std::string>& components, const nlohmann::json& flags) {
+ eos.l_i = eos.l.cast<int>();
- // Allocate the matrix with default models
- std::vector<std::vector<MultiFluidDepartureFunction>> funcs(components.size()); for (auto i = 0; i < funcs.size(); ++i) { funcs[i].resize(funcs.size()); }
+ if (((eos.l_i.cast<double>() - eos.l).cwiseAbs() > 0.0).any()) {
+ throw std::invalid_argument("Non-integer entry in l found");
+ }
- auto depcollection = nlohmann::json::parse(std::ifstream(coolprop_root + "/dev/mixtures/mixture_departure_functions.json"));
+ return eos;
+ };
- auto get_departure_function = [&depcollection](const std::string& Name) {
- for (auto& el : depcollection) {
- if (el["Name"] == Name) { return el; }
+ auto build_gaussian = [&](auto &term) {
+ GaussianEOSTerm eos;
+ eos.n = toeig(term["n"]);
+ eos.t = toeig(term["t"]);
+ eos.d = toeig(term["d"]);
+ eos.eta = toeig(term["eta"]);
+ eos.beta = toeig(term["beta"]);
+ eos.gamma = toeig(term["gamma"]);
+ eos.epsilon = toeig(term["epsilon"]);
+ if (!all_same_length(term, { "n","t","d","eta","beta","gamma","epsilon" })) {
+ throw std::invalid_argument("Lengths are not all identical in Gaussian term");
}
- throw std::invalid_argument("Bad argument");
+ return eos;
+ };
+ auto build_GERG2004 = [&](const auto& term, auto &dep) {
+ if (!all_same_length(term, { "n","t","d","eta","beta","gamma","epsilon" })) {
+ throw std::invalid_argument("Lengths are not all identical in Gaussian term");
+ }
+ auto Npower = term["Npower"];
+ auto NGERG = term["n"].size()- Npower;
+
+ PowerEOSTerm eos;
+ eos.n = toeig(term["n"]).head(Npower);
+ eos.t = toeig(term["t"]).head(Npower);
+ eos.d = toeig(term["d"]).head(Npower);
+ if (term.contains("l")) {
+ eos.l = toeig(term["l"]).head(Npower);
+ }
+ else {
+ eos.l = 0.0 * eos.n;
+ }
+ eos.c = (eos.l > 0).cast<int>().cast<double>();
+ eos.l_i = eos.l.cast<int>();
+ dep.add_term(eos);
+
+ GERG2004EOSTerm e;
+ e.n = toeig(term["n"]).tail(NGERG);
+ e.t = toeig(term["t"]).tail(NGERG);
+ e.d = toeig(term["d"]).tail(NGERG);
+ e.eta = toeig(term["eta"]).tail(NGERG);
+ e.beta = toeig(term["beta"]).tail(NGERG);
+ e.gamma = toeig(term["gamma"]).tail(NGERG);
+ e.epsilon = toeig(term["epsilon"]).tail(NGERG);
+ dep.add_term(e);
+ };
+ auto get_function = [&](auto& funcname) {
+ auto j = get_departure_json(funcname);
+ auto type = j["type"];
+ DepartureTerms dep;
+ if (type == "Exponential") {
+ dep.add_term(build_power(j));
+ }
+ else if(type == "GERG-2004" || type == "GERG-2008") {
+ build_GERG2004(j, dep);
+ }
+ else {
+ throw std::invalid_argument("Bad term type, should not get here");
+ }
+ return dep;
};
for (auto i = 0; i < funcs.size(); ++i) {
for (auto j = i + 1; j < funcs.size(); ++j) {
auto BIP = MultiFluidReducingFunction::get_BIPdep(BIPcollection, { components[i], components[j] }, flags);
- auto function = BIP["function"];
- if (!function.empty()) {
-
- auto info = get_departure_function(function);
- auto N = info["n"].size();
-
- auto toeig = [](const std::vector<double>& v) -> Eigen::ArrayXd { return Eigen::Map<const Eigen::ArrayXd>(&(v[0]), v.size()); };
- auto eigorempty = [&info, &toeig, &N](const std::string& name) -> Eigen::ArrayXd {
- if (!info[name].empty()) {
- return toeig(info[name]);
- }
- else {
- return Eigen::ArrayXd::Zero(N);
- }
- };
-
- MultiFluidDepartureFunction f;
- f.set_type(info["type"]);
- f.n = toeig(info["n"]);
- f.t = toeig(info["t"]);
- f.d = toeig(info["d"]);
-
- f.eta = eigorempty("eta");
- f.beta = eigorempty("beta");
- f.gamma = eigorempty("gamma");
- f.epsilon = eigorempty("epsilon");
-
- Eigen::ArrayXd c(f.n.size()), l(f.n.size()); c.setZero();
- if (info["l"].empty()) {
- // exponential part not included
- l.setZero();
- }
- else {
- l = toeig(info["l"]);
- // l is included, use it to build c; c_i = 1 if l_i > 0, zero otherwise
- for (auto i = 0; i < c.size(); ++i) {
- if (l[i] > 0) {
- c[i] = 1.0;
- }
- }
- }
-
- f.l = l;
- f.c = c;
- funcs[i][j] = f;
- funcs[j][i] = f;
- int rr = 0;
+ std::string funcname = BIP["function"];
+ if (!funcname.empty()) {
+ funcs[i][j] = get_function(funcname);
+ funcs[j][i] = get_function(funcname);
}
else {
- MultiFluidDepartureFunction f;
- f.set_type("none");
- funcs[i][j] = f;
- funcs[j][i] = f;
+ funcs[i][j].add_term(NullEOSTerm());
+ funcs[j][i].add_term(NullEOSTerm());
}
}
}
@@ -447,25 +469,6 @@ auto pow(const mcx::MultiComplex<T> &x, const Eigen::ArrayXd& e) {
return o;
}
-template<class T>
-struct PowIUnaryFunctor {
- const T m_base;
- PowIUnaryFunctor(T base) : m_base(base) {};
- typedef T result_type;
- result_type operator()(const int& e) const{
- switch (e) {
- case 0:
- return 1.0;
- case 1:
- return m_base;
- case 2:
- return m_base * m_base;
- default:
- return powi(m_base, e);
- }
- }
-};
-
inline auto get_EOS_terms(const std::string& coolprop_root, const std::string& name)
{
using namespace nlohmann;
@@ -486,14 +489,6 @@ inline auto get_EOS_terms(const std::string& coolprop_root, const std::string& n
}
}
- auto toeig = [](const std::vector<double>& v) -> Eigen::ArrayXd { return Eigen::Map<const Eigen::ArrayXd>(&(v[0]), v.size()); };
-
- auto all_same_length = [](const nlohmann::json& j, const std::vector<std::string>& ks) {
- std::set<int> lengths;
- for (auto k : ks) { lengths.insert(j[k].size()); }
- return lengths.size() == 1;
- };
-
EOSTerms container;
auto build_power = [&](auto term) {
@@ -501,7 +496,7 @@ inline auto get_EOS_terms(const std::string& coolprop_root, const std::string& n
PowerEOSTerm eos;
- auto eigorzero = [&term, &toeig, &N](const std::string& name) -> Eigen::ArrayXd {
+ auto eigorzero = [&term, &N](const std::string& name) -> Eigen::ArrayXd {
if (!term[name].empty()) {
return toeig(term[name]);
}
diff --git a/include/teqp/models/multifluid_eosterms.hpp b/include/teqp/models/multifluid_eosterms.hpp
index 80547d2..0b6c08b 100644
--- a/include/teqp/models/multifluid_eosterms.hpp
+++ b/include/teqp/models/multifluid_eosterms.hpp
@@ -38,9 +38,22 @@ public:
}
};
+/**
+\f$ \alpha^r = \displaystyle\sum_i n_i \tau^{t_i}\delta^ {d_i} \exp(-\eta_i(\delta-\epsilon_i)^2 -\beta_i(\delta-\gamma_i) }\f$
+*/
+class GERG2004EOSTerm {
+public:
+ Eigen::ArrayXd n, t, d, eta, beta, gamma, epsilon;
+
+ template<typename TauType, typename DeltaType>
+ auto alphar(const TauType& tau, const DeltaType& delta) const {
+ return forceeval((n * exp(t * log(tau) + d * log(delta) - eta * (delta - epsilon).square() - beta * (tau - gamma).square())).sum());
+ }
+};
+
/**
-\f$ \alpha^ r = \displaystyle\sum_i n_i \delta^ { d_i } \tau^ { t_i } \exp(-\delta^ { l_i } - \tau^ { m_i })\f$
+\f$ \alpha^r = \displaystyle\sum_i n_i \delta^ { d_i } \tau^ { t_i } \exp(-\delta^ { l_i } - \tau^ { m_i })\f$
*/
class Lemmon2005EOSTerm {
public:
@@ -54,7 +67,7 @@ public:
};
/**
-\f$ \alpha^ r = \displaystyle\sum_i n_i \tau^{t_i}\delta^ {d_i} \exp(-\eta_i(\delta-\epsilon_i)^2 + \frac{1}{\beta_i(\tau-\gamma_i)^2+b_i}\f$
+\f$ \alpha^r = \displaystyle\sum_i n_i \tau^{t_i}\delta^ {d_i} \exp(-\eta_i(\delta-\epsilon_i)^2 + \frac{1}{\beta_i(\tau-\gamma_i)^2+b_i}\f$
*/
class GaoBEOSTerm {
public:
@@ -67,6 +80,17 @@ public:
}
};
+/**
+\f$ \alpha^r = 0\f$
+*/
+class NullEOSTerm {
+public:
+ template<typename TauType, typename DeltaType>
+ auto alphar(const TauType& tau, const DeltaType& delta) const {
+ return static_cast<std::common_type_t<TauType, DeltaType>>(0.0);
+ }
+};
+
class NonAnalyticEOSTerm {
public:
Eigen::ArrayXd A, B, C, D, a, b, beta, n;
@@ -81,17 +105,24 @@ public:
auto theta = ((1.0 - tau) + A * pow(delta_min1_sq, k)).eval();
auto Delta = (theta.square() + B * pow(delta_min1_sq, a)).eval();
- return forceeval((n * pow(Delta, b) * delta * Psi).eval().sum());
+ auto outval = forceeval((n * pow(Delta, b) * delta * Psi).eval().sum());
+
+ // If we are really, really close to the critical point (tau=delta=1), then the term will become undefined, so let's just return 0 in that case
+ double dbl = getbaseval(outval);
+ if (std::isfinite(dbl)) {
+ return outval;
+ }
+ else {
+ return static_cast<decltype(outval)>(0.0);
+ }
}
};
-
template<typename... Args>
-class EOSTermContainer {
-public:
- using varEOSTerms = std::variant<Args...>;
+class EOSTermContainer {
private:
+ using varEOSTerms = std::variant<Args...>;
std::vector<varEOSTerms> coll;
public:
@@ -107,18 +138,12 @@ public:
std::common_type_t <Tau, Delta> ar = 0.0;
for (const auto& term : coll) {
auto contrib = std::visit([&](auto& t) { return t.alphar(tau, delta); }, term);
- if (std::holds_alternative<NonAnalyticEOSTerm>(term)) {
- double dbl = getbaseval(abs(contrib));
- if (std::isfinite(dbl)) {
- ar = ar + contrib;
- }
- }
- else {
- ar = ar + contrib;
- }
+ ar = ar + contrib;
}
return ar;
}
};
-using EOSTerms = EOSTermContainer<PowerEOSTerm, GaussianEOSTerm, NonAnalyticEOSTerm, Lemmon2005EOSTerm, GaoBEOSTerm, ExponentialEOSTerm>;
\ No newline at end of file
+using EOSTerms = EOSTermContainer<PowerEOSTerm, GaussianEOSTerm, NonAnalyticEOSTerm, Lemmon2005EOSTerm, GaoBEOSTerm, ExponentialEOSTerm>;
+
+using DepartureTerms = EOSTermContainer<PowerEOSTerm, GaussianEOSTerm, GERG2004EOSTerm, NullEOSTerm>;
\ No newline at end of file
diff --git a/include/teqp/types.hpp b/include/teqp/types.hpp
index f315d21..a9644a9 100644
--- a/include/teqp/types.hpp
+++ b/include/teqp/types.hpp
@@ -1,7 +1,10 @@
#pragma once
+#include "nlohmann/json.hpp"
+
#include <vector>
#include <valarray>
+#include <set>
// Registration of types that are considered to be containers
// See https://stackoverflow.com/a/12045843
@@ -28,14 +31,29 @@ auto forceeval(T&& expr)
}
}
+// See https://stackoverflow.com/a/41438758
+template<typename T> struct is_complex_t : public std::false_type {};
+template<typename T> struct is_complex_t<std::complex<T>> : public std::true_type {};
+
template<typename T>
-auto getbaseval(T&& expr)
+auto getbaseval(const T& expr)
{
using namespace autodiff::detail;
if constexpr (isDual<T> || isExpr<T>) {
return val(expr);
}
+ else if constexpr (is_complex_t<T>()) {
+ return expr.real();
+ }
else {
return expr;
}
-}
\ No newline at end of file
+}
+
+auto toeig = [](const std::vector<double>& v) -> Eigen::ArrayXd { return Eigen::Map<const Eigen::ArrayXd>(&(v[0]), v.size()); };
+
+auto all_same_length = [](const nlohmann::json& j, const std::vector<std::string>& ks) {
+ std::set<int> lengths;
+ for (auto k : ks) { lengths.insert(j[k].size()); }
+ return lengths.size() == 1;
+};
--
GitLab