diff --git a/externals/mcx b/externals/mcx
index 9f890c0d7c6a0add8a9afac987e7b3e0eef2c24f..11d97ed4702e144cd9cf125837d8f9d092a984df 160000
--- a/externals/mcx
+++ b/externals/mcx
@@ -1 +1 @@
-Subproject commit 9f890c0d7c6a0add8a9afac987e7b3e0eef2c24f
+Subproject commit 11d97ed4702e144cd9cf125837d8f9d092a984df
diff --git a/include/teqp/critical_tracing.hpp b/include/teqp/critical_tracing.hpp
index a618f2c07d16e44687f986302b9feb99e9ed44b8..6ffd4ffd22e927557fabd2a1a3e741bea9f6b05a 100644
--- a/include/teqp/critical_tracing.hpp
+++ b/include/teqp/critical_tracing.hpp
@@ -62,8 +62,7 @@ auto eigen_problem(const Model& model, const TType T, const RhoType& rhovec) {
}
}
Eigen::MatrixXd Hprime = H(indicesToKeep, indicesToKeep);
- std::cout << H << std::endl;
- std::cout << Hprime << std::endl;
+
auto [eigenvalues, eigenvectors] = sorted_eigen(Hprime);
// Inject values into the U^T and v0 vectors
diff --git a/include/teqp/models/multifluid.hpp b/include/teqp/models/multifluid.hpp
index 8443499399b126b10b2bb1ebbd3c27bc3a72fe0d..444881611b0cbf70046d2c2c8a78d54125277ab8 100644
--- a/include/teqp/models/multifluid.hpp
+++ b/include/teqp/models/multifluid.hpp
@@ -5,6 +5,22 @@
#include <fstream>
#include <string>
+// See https://eigen.tuxfamily.org/dox/TopicCustomizing_CustomScalar.html
+namespace Eigen {
+ template<typename TN> struct NumTraits<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 {
@@ -15,7 +31,7 @@ public:
template<typename TauType, typename DeltaType, typename MoleFractions>
auto alphar(const TauType& tau, const DeltaType& delta, const MoleFractions& molefracs) const {
- using resulttype = decltype(tau* delta* molefracs[0]);
+ using resulttype = std::remove_const<decltype(tau* delta* molefracs[0])>::type; // Type promotion, without the const-ness
resulttype alphar = 0.0;
auto N = molefracs.size();
for (auto i = 0; i < N; ++i) {
@@ -36,12 +52,12 @@ public:
template<typename TauType, typename DeltaType, typename MoleFractions>
auto alphar(const TauType& tau, const DeltaType& delta, const MoleFractions& molefracs) const {
- using resulttype = decltype(tau* delta* molefracs[0]);
+ using resulttype = std::remove_const<decltype(tau* delta* molefracs[0])>::type; // 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);
+ alphar = alphar + molefracs[i] * molefracs[j] * F(i, j) * funcs[i][j].alphar(tau, delta);
}
}
return alphar;
@@ -71,7 +87,7 @@ public:
auto rhored = redfunc.get_rhor(molefrac);
auto delta = rhotot_ / rhored;
auto tau = Tred / T;
- using resulttype = decltype(T* rhovec[0]);
+ using resulttype = decltype(T*rhovec[0]);
return corr.alphar(tau, delta, molefrac) + dep.alphar(tau, delta, molefrac);
}
};
@@ -127,7 +143,7 @@ public:
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*z[i]*z[j]*(z[i] + z[j])/(square(beta(i, j))*z[i] + z[j])*Yij(i, j);
+ sum2 = sum2 + 2.0*z[i]*z[j]*(z[i] + z[j])/(square(beta(i, j))*z[i] + z[j])*Yij(i, j);
}
}
@@ -323,6 +339,24 @@ auto get_departure_function_matrix(const std::string& coolprop_root, const nlohm
return funcs;
}
+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 MultiComplex<T> &x, const Eigen::ArrayXd& e) {
+ Eigen::Array<MultiComplex<T>, Eigen::Dynamic, 1> o(e.size());
+ for (auto i = 0; i < e.size(); ++i) {
+ o[i] = pow(x, e[i]);
+ }
+ return o;
+}
+
class MultiFluidEOS {
public:
enum class types { NOTSETTYPE, GERG2004, GaussianExponential };
diff --git a/src/multifluid.cpp b/src/multifluid.cpp
index c1fc2644010a4365cfc61f048b0fccaa63c77219..3d33fdda55d3bee2d7e80d759c848bec1b45f164 100644
--- a/src/multifluid.cpp
+++ b/src/multifluid.cpp
@@ -22,30 +22,35 @@ auto build_multifluid_model(const std::vector<std::string>& components) {
);
}
-//void trace() {
-// auto model = build_multifluid_model({ "methane", "ethane" });
-// auto rhoc0 = 1.0/model.redfunc.vc[0];
-// auto T = model.redfunc.Tc[0];
-// const auto dT = 1;
-// std::valarray<double> rhovec = { rhoc0, 0.0 };
-// for (auto iter = 0; iter < 1000; ++iter) {
-// auto drhovecdT = get_drhovec_dT_crit(model, T, rhovec);
-// rhovec += drhovecdT * dT;
-// T += dT;
-// int rr = 0;
-// auto z0 = rhovec[0] / rhovec.sum();
-// std::cout << z0 << " ," << rhovec[0] << "," << T << std::endl;
-// if (z0 < 0) {
-// break;
-// }
-// }
-//}
+void trace() {
+ auto model = build_multifluid_model({ "methane", "ethane" });
+ auto rhoc0 = 1.0/model.redfunc.vc[0];
+ auto T = model.redfunc.Tc[0];
+ const auto dT = 1;
+ std::valarray<double> rhovec = { rhoc0, 0.0 };
+ for (auto iter = 0; iter < 1000; ++iter) {
+ auto drhovecdT = get_drhovec_dT_crit(model, T, rhovec);
+ rhovec += drhovecdT * dT;
+ T += dT;
+ int rr = 0;
+ auto z0 = rhovec[0] / rhovec.sum();
+ std::cout << z0 << " ," << rhovec[0] << "," << T << std::endl;
+ if (z0 < 0) {
+ break;
+ }
+ }
+}
int main(){
//test_dummy();
- //trace();
+ trace();
auto model = build_multifluid_model({ "methane", "ethane" });
std::valarray<double> rhovec = { 1.0, 2.0 };
- auto alphar = model.alphar(300.0, rhovec);
+ double T = 300;
+ auto alphar = model.alphar(T, rhovec);
+ double h = 1e-100;
+ auto alpharcom = model.alphar(std::complex<double>(T, h), rhovec).imag()/h;
+ MultiComplex<double> Th{{T, h}};
+ auto alpharcom2 = model.alphar(Th, rhovec).complex().imag()/h;
return EXIT_SUCCESS;
}