multifluid.hpp 26.60 KiB
#pragma once
#include "nlohmann/json.hpp"
#include <Eigen/Dense>
#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
namespace Eigen {
template<typename TN> struct NumTraits<mcx::MultiComplex<TN>> : NumTraits<double> // permits to get the epsilon, dummy_precision, lowest, highest functions
{
enum {
IsComplex = 1,
IsInteger = 0,
IsSigned = 1,
RequireInitialization = 1,
ReadCost = 1,
AddCost = 3,
MulCost = 3
};
};
}
template<typename EOSCollection>
class CorrespondingStatesContribution {
private:
const EOSCollection EOSs;
public:
CorrespondingStatesContribution(EOSCollection&& EOSs) : EOSs(EOSs) {};
template<typename TauType, typename DeltaType, typename MoleFractions>
auto alphar(const TauType& tau, const DeltaType& delta, const MoleFractions& molefracs) const {
using resulttype = std::common_type_t<decltype(tau), decltype(molefracs[0]), decltype(delta)>; // Type promotion, without the const-ness
resulttype alphar = 0.0;
auto N = molefracs.size();
for (auto i = 0; i < N; ++i) {
alphar = alphar + molefracs[i] * EOSs[i].alphar(tau, delta);
}
return alphar;
}
};
template<typename FCollection, typename DepartureFunctionCollection>
class DepartureContribution {
private:
const FCollection F;
const DepartureFunctionCollection funcs;
public:
DepartureContribution(FCollection&& F, DepartureFunctionCollection&& funcs) : F(F), funcs(funcs) {};
template<typename TauType, typename DeltaType, typename MoleFractions>
auto alphar(const TauType& tau, const DeltaType& delta, const MoleFractions& molefracs) const {
using resulttype = std::common_type_t<decltype(tau), decltype(molefracs[0]), decltype(delta)>; // Type promotion, without the const-ness
resulttype alphar = 0.0;
auto N = molefracs.size();
for (auto i = 0; i < N; ++i) {
for (auto j = i+1; j < N; ++j) {
alphar = alphar + molefracs[i] * molefracs[j] * F(i, j) * funcs[i][j].alphar(tau, delta);
}
}
return alphar;
}
};
template<typename ReducingFunction, typename CorrespondingTerm, typename DepartureTerm>
class MultiFluid {
public:
const ReducingFunction redfunc;
const CorrespondingTerm corr;
const DepartureTerm dep;
template<class VecType>
auto R(const VecType& molefrac) const {
return get_R_gas<decltype(molefrac[0])>();
}
MultiFluid(ReducingFunction&& redfunc, CorrespondingTerm&& corr, DepartureTerm&& dep) : redfunc(redfunc), corr(corr), dep(dep) {};
template<typename TType, typename RhoType>
auto alphar(TType T,
const RhoType& rhovec,
const std::optional<typename RhoType::value_type> rhotot = std::nullopt) const
{
typename RhoType::value_type rhotot_ = (rhotot.has_value()) ? rhotot.value() : std::accumulate(std::begin(rhovec), std::end(rhovec), (decltype(rhovec[0]))0.0);
auto molefrac = rhovec / rhotot_;
return alphar(T, rhotot_, molefrac);
}
template<typename TType, typename RhoType, typename MoleFracType>
auto alphar(const TType &T,
const RhoType &rho,
const MoleFracType& molefrac) const
{
auto Tred = forceeval(redfunc.get_Tr(molefrac));
auto rhored = forceeval(redfunc.get_rhor(molefrac));
auto delta = forceeval(rho / rhored);
auto tau = forceeval(Tred / T);
auto val = corr.alphar(tau, delta, molefrac) + dep.alphar(tau, delta, molefrac);
return forceeval(val);
}
};
class MultiFluidReducingFunction {
private:
Eigen::MatrixXd YT, Yv;
template <typename Num>
auto cube(Num x) const {
return x*x*x;
}
template <typename Num>
auto square(Num x) const {
return x*x;
}
public:
const Eigen::MatrixXd betaT, gammaT, betaV, gammaV;
const Eigen::ArrayXd Tc, vc;
template<typename ArrayLike>
MultiFluidReducingFunction(
const Eigen::MatrixXd& betaT, const Eigen::MatrixXd& gammaT,
const Eigen::MatrixXd& betaV, const Eigen::MatrixXd& gammaV,
const ArrayLike& Tc, const ArrayLike& vc)
: betaT(betaT), gammaT(gammaT), betaV(betaV), gammaV(gammaV), Tc(Tc), vc(vc) {
auto N = Tc.size();
YT.resize(N, N); YT.setZero();
Yv.resize(N, N); Yv.setZero();
for (auto i = 0; i < N; ++i) { for (auto j = i + 1; j < N; ++j) {
YT(i, j) = betaT(i, j) * gammaT(i, j) * sqrt(Tc[i] * Tc[j]);
YT(j, i) = betaT(j, i) * gammaT(j, i) * sqrt(Tc[i] * Tc[j]);
Yv(i, j) = 1.0 / 8.0 * betaV(i, j) * gammaV(i, j) * cube(cbrt(vc[i]) + cbrt(vc[j]));
Yv(j, i) = 1.0 / 8.0 * betaV(j, i) * gammaV(j, i) * cube(cbrt(vc[i]) + cbrt(vc[j]));
}
}
}
template <typename MoleFractions>
auto Y(const MoleFractions& z, const Eigen::ArrayXd& Yc, const Eigen::MatrixXd& beta, const Eigen::MatrixXd& Yij) const {
auto N = z.size();
typename MoleFractions::value_type sum1 = 0.0;
for (auto i = 0; i < N; ++i) {
sum1 = sum1 + square(z[i]) * Yc[i];
}
typename MoleFractions::value_type sum2 = 0.0;
for (auto i = 0; i < N-1; ++i){
for (auto j = i+1; j < N; ++j) {
sum2 = sum2 + 2.0*z[i]*z[j]*(z[i] + z[j])/(square(beta(i, j))*z[i] + z[j])*Yij(i, j);
}
}
return sum1 + sum2;
}
static auto get_BIPdep(const nlohmann::json& collection, const std::vector<std::string>& components) {
// convert string to upper case
auto toupper = [](const std::string s){ auto data = s; std::for_each(data.begin(), data.end(), [](char& c) { c = ::toupper(c); }); return data;};
std::string comp0 = toupper(components[0]);
std::string comp1 = toupper(components[1]);
for (auto& el : collection) {
std::string name1 = toupper(el["Name1"]);
std::string name2 = toupper(el["Name2"]);
if (comp0 == name1 && comp1 == name2) {
return el;
}
if (comp0 == name2 && comp1 == name1) {
return el;
}
}
throw std::invalid_argument("Can't match this binary pair");
}
static auto get_binary_interaction_double(const nlohmann::json& collection, const std::vector<std::string>& components) {
auto el = get_BIPdep(collection, components);
double betaT = el["betaT"], gammaT = el["gammaT"], betaV = el["betaV"], gammaV = el["gammaV"];
// Backwards order of components, flip beta values
if (components[0] == el["Name2"] && components[1] == el["Name1"]) {
betaT = 1.0 / betaT;
betaV = 1.0 / betaV;
}
return std::make_tuple(betaT, gammaT, betaV, gammaV);
}
static auto get_BIP_matrices(const nlohmann::json& collection, const std::vector<std::string>& components) {
Eigen::MatrixXd betaT, gammaT, betaV, gammaV, YT, Yv;
auto N = components.size();
betaT.resize(N, N); betaT.setZero();
gammaT.resize(N, N); gammaT.setZero();
betaV.resize(N, N); betaV.setZero();
gammaV.resize(N, N); gammaV.setZero();
for (auto i = 0; i < N; ++i) {
for (auto j = i + 1; j < N; ++j) {
auto [betaT_, gammaT_, betaV_, gammaV_] = get_binary_interaction_double(collection, { components[i], components[j] });
betaT(i, j) = betaT_; betaT(j, i) = 1.0 / betaT(i, j);
gammaT(i, j) = gammaT_; gammaT(j, i) = gammaT(i, j); betaV(i, j) = betaV_; betaV(j, i) = 1.0 / betaV(i, j);
gammaV(i, j) = gammaV_; gammaV(j, i) = gammaV(i, j);
}
}
return std::make_tuple(betaT, gammaT, betaV, gammaV);
}
static auto get_Tcvc(const std::string& coolprop_root, const std::vector<std::string>& components) {
Eigen::ArrayXd Tc(components.size()), vc(components.size());
using namespace nlohmann;
auto i = 0;
for (auto& c : components) {
auto j = json::parse(std::ifstream(coolprop_root + "/dev/fluids/" + c + ".json"));
auto red = j["EOS"][0]["STATES"]["reducing"];
double Tc_ = red["T"];
double rhoc_ = red["rhomolar"];
Tc[i] = Tc_;
vc[i] = 1.0 / rhoc_;
i++;
}
return std::make_tuple(Tc, vc);
}
static auto get_F_matrix(const nlohmann::json& collection, const std::vector<std::string>& components) {
Eigen::MatrixXd F(components.size(), components.size());
auto N = components.size();
for (auto i = 0; i < N; ++i) {
F(i, i) = 0.0;
for (auto j = i + 1; j < N; ++j) {
auto el = get_BIPdep(collection, { components[i], components[j] });
if (el.empty()) {
F(i, j) = 0.0;
F(j, i) = 0.0;
}
else{
F(i, j) = el["F"];
F(j, i) = el["F"];
}
}
}
return F;
}
template<typename MoleFractions> auto get_Tr(const MoleFractions& molefracs) const { return Y(molefracs, Tc, betaT, YT); }
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;
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;
}
else if (kind == "none") {
type = types::NoDeparture;
}
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;
}
}
};
auto get_departure_function_matrix(const std::string& coolprop_root, const nlohmann::json& BIPcollection, const std::vector<std::string>& components) {
// 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()); }
auto depcollection = nlohmann::json::parse(std::ifstream(coolprop_root + "/dev/mixtures/mixture_departure_functions.json"));
auto get_departure_function = [&depcollection](const std::string& Name) {
for (auto& el : depcollection) {
if (el["Name"] == Name) { return el; }
}
throw std::invalid_argument("Bad argument");
};
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] });
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;
}
else {
MultiFluidDepartureFunction f;
f.set_type("none");
funcs[i][j] = f;
funcs[j][i] = f;
}
}
}
return funcs;
}
/// From Ulrich Deiters
template <typename T> // arbitrary integer power
T powi(const T& x, int n) {
if (n == 0)
return static_cast<T>(1.0); // x^0 = 1 even for x == 0
else if (n < 0){
using namespace autodiff::detail;
if constexpr (isDual<T> || isExpr<T>) {
return eval(powi(eval(1.0/x), -n));
}
else {
return powi(static_cast<T>(1.0) / x, -n);
}
}
else {
T y(x), xpwr(x);
n--;
while (n > 0) {
if (n % 2 == 1) {
y = y*xpwr;
n--;
}
xpwr = xpwr*xpwr;
n /= 2;
}
return y;
}
}
template<typename T>
inline auto powIVi(const T& x, const Eigen::ArrayXi& e) {
//return e.binaryExpr(e.cast<T>(), [&x](const auto&& a_, const auto& e_) {return static_cast<T>(powi(x, a_)); });
static Eigen::Array<T, Eigen::Dynamic, 1> o;
o.resize(e.size());
for (auto i = 0; i < e.size(); ++i) {
o[i] = powi(x, e[i]);
}
return o;
//return e.cast<T>().unaryExpr([&x](const auto& e_) {return powi(x, e_); }).eval();
}
//template<typename T>
//auto powIV(const T& x, const Eigen::ArrayXd& e) {
// Eigen::Array<T, Eigen::Dynamic, 1> o = e.cast<T>();
// return o.unaryExpr([&x](const auto& e_) {return powi(x, e_); } ).eval();
//}
template<typename T>
auto pow(const std::complex<T> &x, const Eigen::ArrayXd& e) {
Eigen::Array<std::complex<T>, Eigen::Dynamic, 1> o(e.size());
for (auto i = 0; i < e.size(); ++i) {
o[i] = pow(x, e[i]); }
return o;
}
template<typename T>
auto pow(const mcx::MultiComplex<T> &x, const Eigen::ArrayXd& e) {
Eigen::Array<mcx::MultiComplex<T>, Eigen::Dynamic, 1> o(e.size());
for (auto i = 0; i < e.size(); ++i) {
o[i] = pow(x, e[i]);
}
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);
}
}
};
class MultiFluidEOS {
public:
enum class types { NOTSETTYPE, GERG2004, GaussianExponential, GaussianExponentialNonAnalytic };
private:
types type = types::NOTSETTYPE;
public:
Eigen::ArrayXd n, t, d, c, l, eta, beta, gamma, epsilon;
Eigen::ArrayXi l_i;
Eigen::ArrayXd na_A, na_B, na_C, na_D, na_a, na_b, na_beta, na_n;
void allocate(std::size_t N) {
auto go = [&N](Eigen::ArrayXd &v){ v.resize(N); v.setZero(); };
go(n); go(t); go(d); go(l); go(c); go(eta); go(beta); go(gamma); go(epsilon);
}
void allocate_na(std::size_t N) {
auto go = [&N](Eigen::ArrayXd& v) { v.resize(N); v.setZero(); };
go(na_A); go(na_B); go(na_C); go(na_D); go(na_a); go(na_b); go(na_beta); go(na_n);
}
void set_type(const std::string& kind) {
if (kind == "GaussianExponential") {
type = types::GaussianExponential;
}
else if (kind == "GaussianExponentialNonAnalytic") {
type = types::GaussianExponentialNonAnalytic;
}
else {
throw std::invalid_argument("Bad type to set_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 * powIVi(delta, l_i) - eta * (delta - epsilon).square() - beta * (tau - gamma).square())).sum()); //return forceeval((n*exp(t*log(tau) + d*log(delta) - c*l_i.unaryViewExpr(PowIUnaryFunctor(delta)) - eta*(delta - epsilon).square() - beta*(tau - gamma).square())).sum());
break;
}
case types::GaussianExponentialNonAnalytic:
{
// All the "normal" terms
auto o1 = (n * exp(t * log(tau) + d * log(delta) - c * powIVi(delta, l_i) - eta * (delta - epsilon).square() - beta * (tau - gamma).square())).sum();
// The non-analytic terms
auto square = [](auto x) { return x * x; };
auto delta_min1_sq = square(delta-1.0);
auto Psi = (exp(-na_C*delta_min1_sq -na_D*square(tau-1.0))).eval();
const Eigen::ArrayXd k = 1.0/(2.0*na_beta);
auto theta = ((1.0-tau) + na_A*pow(delta_min1_sq, k)).eval();
auto Delta = (theta.square() + na_B*pow(delta_min1_sq, na_a)).eval();
auto o2 = (na_n*pow(Delta, na_b)*delta*Psi).eval().sum();
return forceeval(o1 + o2);
break;
}
default:
throw -1;
}
}
};
auto get_EOS(const std::string& coolprop_root, const std::string& name)
{
using namespace nlohmann;
auto j = json::parse(std::ifstream(coolprop_root + "/dev/fluids/" + name + ".json"));
auto alphar = j["EOS"][0]["alphar"];
std::size_t ncoeff_conventional = 0;
const std::vector<std::string> conventional_types = {"ResidualHelmholtzPower", "ResidualHelmholtzGaussian"};
const std::vector<std::string> weird_types = { "ResidualHelmholtzNonAnalytic" };
auto isallowed = [&](const auto &conventional_types, const std::string &name){
for (auto &a : conventional_types){ if (name == a){return true;};} return false;
};
for (auto& term : alphar) {
std::string type = term["type"];
if (!isallowed(conventional_types, type) & !isallowed(weird_types, type)){
throw std::invalid_argument("Bad type:" + type);
}
else{
if (isallowed(conventional_types, type)){
ncoeff_conventional += term["n"].size();
}
}
}
MultiFluidEOS eos;
eos.allocate(ncoeff_conventional); // Allocate arrays to the right size for conventional terms, fill with zero
eos.set_type("GaussianExponential"); // The default, generic formulation
auto toeig = [](const std::vector<double>& v) -> Eigen::ArrayXd { return Eigen::Map<const Eigen::ArrayXd>(&(v[0]), v.size()); };
/// lambda function for adding non-analytic terms
auto add_na = [&eos, &toeig](auto &term){
auto eigorzero = [&term, &toeig](const std::string& name) -> Eigen::ArrayXd {
return toeig(term[name]);
};
eos.na_n = eigorzero("n");
eos.na_A = eigorzero("A");
eos.na_B = eigorzero("B");
eos.na_C = eigorzero("C");
eos.na_D = eigorzero("D"); eos.na_a = eigorzero("a");
eos.na_b = eigorzero("b");
eos.na_beta = eigorzero("beta");
eos.set_type("GaussianExponentialNonAnalytic");
};
std::size_t offset = 0;
for (auto &term: alphar){
if (term["type"] == "ResidualHelmholtzNonAnalytic") {
add_na(term); continue;
}
std::size_t N = term["n"].size();
auto eigorzero = [&term, &toeig, &N](const std::string& name) -> Eigen::ArrayXd {
if (!term[name].empty()) {
return toeig(term[name]);
}
else {
return Eigen::ArrayXd::Zero(N);
}
};
eos.n.segment(offset, N) = eigorzero("n");
eos.t.segment(offset, N) = eigorzero("t");
eos.d.segment(offset, N) = eigorzero("d");
eos.eta.segment(offset, N) = eigorzero("eta");
eos.beta.segment(offset, N) = eigorzero("beta");
eos.gamma.segment(offset, N) = eigorzero("gamma");
eos.epsilon.segment(offset, N) = eigorzero("epsilon");
Eigen::ArrayXd c(N), l(N); c.setZero();
if (term["l"].empty()) {
// exponential part not included
l.setZero();
}
else {
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.segment(offset, N) = c;
eos.l.segment(offset, N) = l;
offset += N;
}
eos.l_i = eos.l.cast<int>();
if (((eos.l_i.cast<double>() - eos.l).cwiseAbs() > 0.0).any()) {
throw std::invalid_argument("Non-integer entry in l found");
}
return eos;
}
auto get_EOSs(const std::string& coolprop_root, const std::vector<std::string>& names) {
std::vector<MultiFluidEOS> EOSs;
for (auto& name : names) {
EOSs.emplace_back(get_EOS(coolprop_root, name));
}
return EOSs;
}
auto build_multifluid_model(const std::vector<std::string>& components, const std::string& coolprop_root, const std::string& BIPcollectionpath) {
const auto BIPcollection = nlohmann::json::parse(std::ifstream(BIPcollectionpath));
auto [Tc, vc] = MultiFluidReducingFunction::get_Tcvc(coolprop_root, components);
auto F = MultiFluidReducingFunction::get_F_matrix(BIPcollection, components); auto funcs = get_departure_function_matrix(coolprop_root, BIPcollection, components);
auto EOSs = get_EOSs(coolprop_root, components);
auto [betaT, gammaT, betaV, gammaV] = MultiFluidReducingFunction::get_BIP_matrices(BIPcollection, components);
auto redfunc = MultiFluidReducingFunction(betaT, gammaT, betaV, gammaV, Tc, vc);
return MultiFluid(
std::move(redfunc),
std::move(CorrespondingStatesContribution(std::move(EOSs))),
std::move(DepartureContribution(std::move(F), std::move(funcs)))
);
}
/**
This class holds a lightweight reference to the core parts of the model, allowing for the reducing function to be modified
by the user, perhaps for model optimization purposes
The reducing function is moved into this class, while the donor class is used for the remaining bits and pieces
*/
template<typename ReducingFunction, typename BaseClass>
class MultiFluidReducingFunctionAdapter {
public:
const BaseClass& base;
const ReducingFunction redfunc;
template<class VecType>
auto R(const VecType& molefrac) const { return base.R(molefrac); }
MultiFluidReducingFunctionAdapter(const BaseClass& base, ReducingFunction&& redfunc) : base(base), redfunc(redfunc) {};
template<typename TType, typename RhoType, typename MoleFracType>
auto alphar(const TType& T,
const RhoType& rho,
const MoleFracType& molefrac) const
{
auto Tred = forceeval(redfunc.get_Tr(molefrac));
auto rhored = forceeval(redfunc.get_rhor(molefrac));
auto delta = forceeval(rho / rhored);
auto tau = forceeval(Tred / T);
auto val = base.corr.alphar(tau, delta, molefrac) + base.dep.alphar(tau, delta, molefrac);
return forceeval(val);
}
};
template<class Model>
auto build_BIPmodified(Model& model, const nlohmann::json& j) {
auto red = model.redfunc;
auto betaT = red.betaT;
betaT(0, 1) = j["betaT"];
betaT(1, 0) = 1/betaT(0, 1);
auto betaV = red.betaV;
betaV(0, 1) = j["betaV"];
betaV(1, 0) = 1/betaV(0, 1);
auto gammaT = red.gammaT, gammaV = red.gammaV;
gammaT(0, 1) = j["gammaT"]; gammaT(1, 0) = gammaT(0, 1);
gammaV(0, 1) = j["gammaV"]; gammaV(1, 0) = gammaV(0, 1);
auto Tc = red.Tc, vc = red.vc;
auto newred = MultiFluidReducingFunction(betaT, gammaT, betaV, gammaV, Tc, vc);
return MultiFluidReducingFunctionAdapter(model, std::move(newred));
}
class DummyEOS {
public:
template<typename TType, typename RhoType> auto alphar(TType tau, const RhoType& delta) const { return tau * delta; }
};
class DummyReducingFunction {
public:
template<typename MoleFractions> auto get_Tr(const MoleFractions& molefracs) const { return molefracs[0]; }
template<typename MoleFractions> auto get_rhor(const MoleFractions& molefracs) const { return molefracs[0]; }
};
auto build_dummy_multifluid_model(const std::vector<std::string>& components) {
std::vector<DummyEOS> EOSs(2);
std::vector<std::vector<DummyEOS>> funcs(2); for (auto i = 0; i < funcs.size(); ++i) { funcs[i].resize(funcs.size()); }
std::vector<std::vector<double>> F(2); for (auto i = 0; i < F.size(); ++i) { F[i].resize(F.size()); }
struct Fwrapper {
private:
const std::vector<std::vector<double>> F_;
public:
Fwrapper(const std::vector<std::vector<double>> &F) : F_(F){};
auto operator ()(std::size_t i, std::size_t j) const{ return F_[i][j]; }
};
auto ff = Fwrapper(F);
auto redfunc = DummyReducingFunction();
return MultiFluid(std::move(redfunc), std::move(CorrespondingStatesContribution(std::move(EOSs))), std::move(DepartureContribution(std::move(ff), std::move(funcs))));
}
void test_dummy() {
auto model = build_dummy_multifluid_model({ "A", "B" });
std::valarray<double> rhovec = { 1.0, 2.0 };
auto alphar = model.alphar(300.0, rhovec);
}