diff --git a/externals/mcx b/externals/mcx
index 5e48688a952d1e24e876541a68609e43a7083aa3..fa339e9963f58fb4998e56271173c9afa3e5dc28 160000
--- a/externals/mcx
+++ b/externals/mcx
@@ -1 +1 @@
-Subproject commit 5e48688a952d1e24e876541a68609e43a7083aa3
+Subproject commit fa339e9963f58fb4998e56271173c9afa3e5dc28
diff --git a/include/teqp/algorithms/VLE.hpp b/include/teqp/algorithms/VLE.hpp
index 70050a45717ae8bfc99504dcb7271101fb2e3292..f9ef0a9d0bbf82c856fc6d2bfa8a2aea4e3ce783 100644
--- a/include/teqp/algorithms/VLE.hpp
+++ b/include/teqp/algorithms/VLE.hpp
@@ -22,11 +22,12 @@ auto linsolve(const A& a, const B& b) {
return a.matrix().colPivHouseholderQr().solve(b.matrix()).array().eval();
}
-template<typename Model, typename TYPE = double>
+template<typename Model, typename TYPE = double, ADBackends backend = ADBackends::autodiff>
class IsothermPureVLEResiduals {
- typedef Eigen::Array<TYPE, 2, 1> EigenArray;
- typedef Eigen::Array<TYPE, 1, 1> EigenArray1;
- typedef Eigen::Array<TYPE, 2, 2> EigenMatrix;
+public:
+ using EigenArray = Eigen::Array<TYPE, 2, 1>;
+ using EigenArray1 = Eigen::Array<TYPE, 1, 1>;
+ using EigenMatrix = Eigen::Array<TYPE, 2, 2>;
private:
const Model& m_model;
TYPE m_T;
@@ -59,13 +60,13 @@ public:
const TYPE R = m_model.R(molefracs);
double R0_over_Rr = R0 / Rr;
- auto derL = tdx::template get_Ar0n<2>(m_model, T, rhomolarL, molefracs);
+ auto derL = tdx::template get_Ar0n<2, backend>(m_model, T, rhomolarL, molefracs);
auto pRTL = rhomolarL*(R0_over_Rr + derL[1]); // p/(R*T)
auto dpRTLdrhoL = R0_over_Rr + 2*derL[1] + derL[2];
auto hatmurL = derL[1] + derL[0] + R0_over_Rr*log(rhomolarL);
auto dhatmurLdrho = (2*derL[1] + derL[2])/rhomolarL + R0_over_Rr/rhomolarL;
- auto derV = tdx::template get_Ar0n<2>(m_model, T, rhomolarV, molefracs);
+ auto derV = tdx::template get_Ar0n<2, backend>(m_model, T, rhomolarV, molefracs);
auto pRTV = rhomolarV*(R0_over_Rr + derV[1]); // p/(R*T)
auto dpRTVdrhoV = R0_over_Rr + 2*derV[1] + derV[2];
auto hatmurV = derV[1] + derV[0] + R0_over_Rr *log(rhomolarV);
@@ -97,8 +98,9 @@ public:
};
template<typename Residual, typename Scalar>
-Eigen::ArrayXd do_pure_VLE_T(Residual &resid, Scalar rhoL, Scalar rhoV, int maxiter) {
- auto rhovec = (Eigen::ArrayXd(2) << rhoL, rhoV).finished();
+auto do_pure_VLE_T(Residual &resid, Scalar rhoL, Scalar rhoV, int maxiter) {
+ using EArray = Eigen::Array<Scalar, 2, 1>;
+ auto rhovec = (EArray() << rhoL, rhoV).finished();
auto r0 = resid.call(rhovec);
auto J = resid.Jacobian(rhovec);
for (int iter = 0; iter < maxiter; ++iter){
@@ -108,21 +110,28 @@ Eigen::ArrayXd do_pure_VLE_T(Residual &resid, Scalar rhoL, Scalar rhoV, int maxi
}
auto v = J.matrix().colPivHouseholderQr().solve(-r0.matrix()).array().eval();
auto rhovecnew = (rhovec + v).eval();
+ double r00 = static_cast<double>(r0[0]);
+ double r01 = static_cast<double>(r0[1]);
// If the solution has stopped improving, stop. The change in rhovec is equal to v in infinite precision, but
// not when finite precision is involved, use the minimum non-denormal float as the determination of whether
// the values are done changing
- if (((rhovecnew - rhovec).cwiseAbs() < std::numeric_limits<Scalar>::min()).all()) {
+ auto minval = std::numeric_limits<Scalar>::epsilon();
+ double minvaldbl = static_cast<double>(minval);
+ if (((rhovecnew - rhovec).cwiseAbs() < minval).all()) {
+ break;
+ }
+ if ((r0.cwiseAbs() < minval).all()) {
break;
}
rhovec = rhovecnew;
}
- return (Eigen::ArrayXd(2) << rhovec[0], rhovec[1]).finished();
+ return rhovec;
}
-template<typename Model, typename Scalar>
-Eigen::ArrayXd pure_VLE_T(const Model& model, Scalar T, Scalar rhoL, Scalar rhoV, int maxiter) {
- auto res = IsothermPureVLEResiduals(model, T);
+template<typename Model, typename Scalar, ADBackends backend = ADBackends::autodiff>
+auto pure_VLE_T(const Model& model, Scalar T, Scalar rhoL, Scalar rhoV, int maxiter) {
+ auto res = IsothermPureVLEResiduals<Model, Scalar, backend>(model, T);
return do_pure_VLE_T(res, rhoL, rhoV, maxiter);
}
diff --git a/include/teqp/derivs.hpp b/include/teqp/derivs.hpp
index e6196337e520aa612a977a8307eb6d26768d5f9e..64b3f6693befa1acd2c029b222c515477c469cbd 100644
--- a/include/teqp/derivs.hpp
+++ b/include/teqp/derivs.hpp
@@ -205,7 +205,7 @@ struct TDXDerivatives {
template<int Nderiv, ADBackends be = ADBackends::autodiff>
static auto get_Ar0n(const Model& model, const Scalar& T, const Scalar& rho, const VectorType& molefrac) {
- std::valarray<double> o(Nderiv+1);
+ std::valarray<Scalar> o(Nderiv+1);
if constexpr (be == ADBackends::autodiff) {
autodiff::HigherOrderDual<Nderiv, double> rhodual = rho;
auto f = [&model, &T, &molefrac](const auto& rho_) { return eval(model.alphar(T, rho_, molefrac)); };
diff --git a/include/teqp/models/cubics.hpp b/include/teqp/models/cubics.hpp
index 76782c0f5fe1328885bfb13f3b599459326d627a..a631d6faf739f6111a10b3f73d5bae82ab2f60ac 100644
--- a/include/teqp/models/cubics.hpp
+++ b/include/teqp/models/cubics.hpp
@@ -141,7 +141,7 @@ public:
throw std::invalid_argument("Sizes do not match");
}
auto b = get_b(T, molefrac);
- auto Psiminus = -log(1.0 - b * rho);
+ auto Psiminus = -1.0*log(1.0 - b * rho);
auto Psiplus = log((Delta1 * b * rho + 1.0) / (Delta2 * b * rho + 1.0)) / (b * (Delta1 - Delta2));
auto val = Psiminus - get_a(T, molefrac) / (Ru * T) * Psiplus;
return forceeval(val);
@@ -178,8 +178,8 @@ auto canonical_SRK(TCType Tc_K, PCType pc_K, AcentricType acentric, const Eigen:
template <typename TCType, typename PCType, typename AcentricType>
auto canonical_PR(TCType Tc_K, PCType pc_K, AcentricType acentric, const Eigen::ArrayXXd& kmat = {}) {
- double Delta1 = 1+sqrt(2);
- double Delta2 = 1-sqrt(2);
+ double Delta1 = 1+sqrt(2.0);
+ double Delta2 = 1-sqrt(2.0);
AcentricType m = acentric*0.0;
std::vector<AlphaFunctionOptions> alphas;
for (auto i = 0; i < Tc_K.size(); ++i) {
diff --git a/src/test_mpderivs.cpp b/src/test_mpderivs.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..80a19a2b0f95fa59cfcb6dc54de13375c344c050
--- /dev/null
+++ b/src/test_mpderivs.cpp
@@ -0,0 +1,27 @@
+// Imports from boost
+#include <boost/multiprecision/cpp_bin_float.hpp>
+using namespace boost::multiprecision;
+
+#include "teqp/models/cubics.hpp"
+#include "teqp/algorithms/VLE.hpp"
+#include "teqp/derivs.hpp"
+
+int main() {
+
+ // Values taken from http://dx.doi.org/10.6028/jres.121.011
+ std::valarray<double> Tc_K = { 154.581}, pc_Pa = { 5042800}, acentric = { 0.022};
+ auto modelPR = teqp::canonical_PR(Tc_K, pc_Pa, acentric);
+
+ using my_float = boost::multiprecision::number<boost::multiprecision::cpp_bin_float<200>>;
+ auto [rhoLdbl, rhoVdbl] = modelPR.superanc_rhoLV(125);
+ my_float T = 125;
+ auto soln = teqp::pure_VLE_T<decltype(modelPR), my_float, teqp::ADBackends::multicomplex>(modelPR, T, static_cast<my_float>(rhoLdbl), static_cast<my_float>(rhoVdbl), 20).cast<double>();
+ std::cout << soln << std::endl;
+
+ my_float x = 10567.6789;
+ std::function<mcx::MultiComplex<my_float>(const mcx::MultiComplex<my_float>&)> ff = [](const auto & x) {return 1.0 - x;};
+ auto o = mcx::diff_mcx1<my_float>(ff, x, 3, true);
+ std::cout << o[1] << std::endl;
+
+ std::cout << std::is_arithmetic_v<my_float> << std::endl;
+}
\ No newline at end of file