From 3779b027c6f2aac6e94ebd19ddc48f4fb9a6364e Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.bell@nist.gov>
Date: Wed, 31 Mar 2021 19:55:46 -0400
Subject: [PATCH] Everything compiles again minus most of the valarray
Isochoric derivatives moved into their own namespace
---
externals/mcx | 2 +-
include/teqp/algorithms/critical_tracing.hpp | 71 +++++-----
include/teqp/derivs.hpp | 138 ++++++++++++-------
src/multifluid.cpp | 61 +++++---
src/tests/catch_tests.cxx | 33 ++++-
5 files changed, 200 insertions(+), 105 deletions(-)
diff --git a/externals/mcx b/externals/mcx
index ffc9c34..1f97ce8 160000
--- a/externals/mcx
+++ b/externals/mcx
@@ -1 +1 @@
-Subproject commit ffc9c34e8ddfcebec9d417096735b4516e673c98
+Subproject commit 1f97ce856ea22bf90c483253fcbc8da2ebc2b9db
diff --git a/include/teqp/algorithms/critical_tracing.hpp b/include/teqp/algorithms/critical_tracing.hpp
index 44baf51..fca4493 100644
--- a/include/teqp/algorithms/critical_tracing.hpp
+++ b/include/teqp/algorithms/critical_tracing.hpp
@@ -32,13 +32,15 @@ auto eigen_problem(const Model& model, const TType T, const RhoType& rhovec) {
EigenData ed;
auto N = rhovec.size();
- std::valarray<bool> mask = (rhovec != 0);
+ Eigen::ArrayX<bool> mask = (rhovec != 0).eval();
+
+ using id = IsochoricDerivatives<decltype(model)>;
// Build the Hessian for the residual part;
#if defined(USE_AUTODIFF)
- auto H = build_Psir_Hessian_autodiff(model, T, rhovec);
+ auto H = id::build_Psir_Hessian_autodiff(model, T, rhovec);
#else
- auto H = build_Psir_Hessian_mcx(model, T, rhovec);
+ auto H = id::build_Psir_Hessian_mcx(model, T, rhovec);
#endif
// ... and add ideal-gas terms to H
for (auto i = 0; i < N; ++i) {
@@ -100,7 +102,7 @@ auto eigen_problem(const Model& model, const TType T, const RhoType& rhovec) {
}
struct psi1derivs {
- std::valarray<double> psir, psi0, tot;
+ Eigen::ArrayXd psir, psi0, tot;
EigenData ei;
};
@@ -112,7 +114,7 @@ auto get_derivs(const Model& model, const TType T, const RhoType& rhovec) {
auto ei = eigen_problem(model, T, rhovec);
// Ideal-gas contributions of psi0 w.r.t. sigma_1, in the same form as the residual part
- std::valarray<double> psi0_derivs(0.0, 5);
+ Eigen::ArrayXd psi0_derivs(5); psi0_derivs.setZero();
psi0_derivs[0] = -1; // Placeholder, not needed
psi0_derivs[1] = -1; // Placeholder, not needed
for (auto i = 0; i < rhovec.size(); ++i) {
@@ -129,23 +131,28 @@ auto get_derivs(const Model& model, const TType T, const RhoType& rhovec) {
VectorXdual4th rhovecad(rhovec.size()); for (auto i = 0; i < rhovec.size(); ++i) { rhovecad[i] = rhovec[i]; }
dual4th varsigma{0.0};
auto wrapper = [&rhovecad, &v0, &T, &model](const auto &sigma_1) {
- auto rhovecused = rhovecad + sigma_1*v0;
- return get_Psir(model, T, rhovecused);
+ auto rhovecused = (rhovecad + sigma_1*v0).eval();
+ auto rhotot = rhovecused.sum();
+ auto molefrac = (rhovecused / rhotot).eval();
+ return eval(model.alphar(T, rhotot, molefrac) * model.R * T * rhotot);
};
- auto derivs = derivatives(wrapper, wrt(varsigma, varsigma, varsigma, varsigma), at(varsigma));
- auto psir_derivs = std::valarray<double>(&derivs[0], derivs.size());
+ auto psir_derivs_ = derivatives(wrapper, wrt(varsigma), at(varsigma));
+ RhoType psir_derivs; psir_derivs.resize(5);
+ for (auto i = 0; i < 5; ++i){ psir_derivs[i] = psir_derivs_[i]; }
+
#else
using namespace mcx;
// Calculate the first through fourth derivative of Psi^r w.r.t. sigma_1
- std::valarray<MultiComplex<double>> v0(ei.v0.size()); for (auto i = 0; i < ei.v0.size(); ++i) { v0[i] = ei.v0[i]; }
- std::valarray<MultiComplex<double>> rhovecmcx(rhovec.size()); for (auto i = 0; i < rhovec.size(); ++i) { rhovecmcx[i] = rhovec[i]; }
+ Eigen::Vector<MultiComplex<double>,Eigen::Dynamic> v0(ei.v0.size()); for (auto i = 0; i < ei.v0.size(); ++i) { v0[i] = ei.v0[i]; }
+ Eigen::Vector<MultiComplex<double>,Eigen::Dynamic> rhovecmcx(rhovec.size()); for (auto i = 0; i < rhovec.size(); ++i) { rhovecmcx[i] = rhovec[i]; }
using fcn_t = std::function<MultiComplex<double>(const MultiComplex<double>&)>;
fcn_t wrapper = [&rhovecmcx, &v0, &T, &model](const MultiComplex<double>& sigma_1) {
- std::valarray<MultiComplex<double>> rhovecused = rhovecmcx + sigma_1 * v0;
- return get_Psir(model, T, rhovecused);
+ Eigen::Vector<MultiComplex<double>, Eigen::Dynamic> rhovecused = rhovecmcx + sigma_1 * v0;
+ auto rhotot = rhovecused.sum();
+ auto molefrac = rhovecused/rhotot;
+ return model.alphar(T, rhotot, molefrac)*model.R*T*rhotot;
};
- auto psir_derivs_ = diff_mcx1(wrapper, 0.0, 4, true);
- auto psir_derivs = std::valarray<double>(&psir_derivs_[0], psir_derivs_.size());
+ RhoType psir_derivs = diff_mcx1(wrapper, 0.0, 4, true);
#endif
// As a sanity check, the minimum eigenvalue of the Hessian constructed based on the molar concentrations
@@ -170,7 +177,7 @@ bool all(const Iterable& foo) {
template <typename Iterable>
bool any(const Iterable& foo) {
return std::any_of(std::begin(foo), std::end(foo), [](const auto x) { return x; });
-}
+}
template<typename Model, typename TType, typename RhoType>
auto get_drhovec_dT_crit(const Model& model, const TType T, const RhoType& rhovec) {
@@ -191,12 +198,11 @@ auto get_drhovec_dT_crit(const Model& model, const TType T, const RhoType& rhove
auto sigma2 = 2e-5 * rhovec.sum(); // This is the perturbation along the second eigenvector
- auto v1 = std::valarray<double>(&ei.v1[0], ei.v1.size());
- decltype(v1) rhovec_plus = rhovec + v1 * sigma2;
- decltype(v1) rhovec_minus = rhovec - v1 * sigma2;
+ auto rhovec_plus = (rhovec + ei.v1 * sigma2).eval();
+ auto rhovec_minus = (rhovec - ei.v1 * sigma2).eval();
std::string stepping_desc = "";
auto deriv_sigma2 = all_derivs.tot;
- auto eval = [](const auto &ex){ return std::valarray<bool>(ex); };
+ auto eval = [](const auto &ex){ return ex.eval(); };
if (all(eval(rhovec_minus > 0)) && all(eval(rhovec_plus > 0))) {
// Conventional centered derivative
auto plus_sigma2 = get_derivs(model, T, rhovec_plus);
@@ -207,7 +213,7 @@ auto get_drhovec_dT_crit(const Model& model, const TType T, const RhoType& rhove
else if (any(eval(rhovec_minus < 0))) {
// Forward derivative in the direction of v1
auto plus_sigma2 = get_derivs(model, T, rhovec_plus);
- auto rhovec_2plus = rhovec + 2 * v1 * sigma2;
+ auto rhovec_2plus = (rhovec + 2 * ei.v1 * sigma2).eval();
auto plus2_sigma2 = get_derivs(model, T, rhovec_2plus);
deriv_sigma2 = (-3 * derivs + 4 * plus_sigma2.tot - plus2_sigma2.tot) / (2.0 * sigma2);
stepping_desc = "forward";
@@ -215,7 +221,7 @@ auto get_drhovec_dT_crit(const Model& model, const TType T, const RhoType& rhove
else if (any(eval(rhovec_minus > 0))) {
// Negative derivative in the direction of v1
auto minus_sigma2 = get_derivs(model, T, rhovec_minus);
- auto rhovec_2minus = rhovec - 2 * v1 * sigma2;
+ auto rhovec_2minus = (rhovec - 2 * ei.v1 * sigma2).eval();
auto minus2_sigma2 = get_derivs(model, T, rhovec_2minus);
deriv_sigma2 = (-3 * derivs + 4 * minus_sigma2.tot - minus2_sigma2.tot) / (-2.0 * sigma2);
stepping_desc = "backwards";
@@ -245,14 +251,14 @@ auto get_drhovec_dT_crit(const Model& model, const TType T, const RhoType& rhove
// print("RHS", RHS)
// print("drhovec_dT", drhovec_dT)
- return std::valarray<double>(&drhovec_dT(0), drhovec_dT.size());
+ return drhovec_dT;
}
template<typename ModelType, typename VecType>
auto critical_polish_molefrac(const ModelType &model, const double T, const VecType &rhovec, const double z0) {
auto polish_x_resid = [&model, &z0](const auto& x) {
auto T = x[0];
- std::valarray<double> rhovec = { x[1], x[2] };
+ Eigen::ArrayXd rhovec(2); rhovec << x[1], x[2] ;
auto z0new = rhovec[0] / rhovec.sum();
auto derivs = get_derivs(model, T, rhovec);
// First two are residuals on critical point, third is residual on composition
@@ -266,14 +272,15 @@ auto critical_polish_molefrac(const ModelType &model, const double T, const VecT
if (!std::isfinite(T) || !std::isfinite(x[1]) || !std::isfinite(x[2])) {
throw std::invalid_argument("Something not finite; aborting polishing");
}
- return std::make_tuple(x[0], x.tail(x.size()-1).eval());
+ Eigen::ArrayXd rho = x.tail(x.size() - 1);
+ return std::make_tuple(x[0], rho);
}
template<typename ModelType, typename VecType>
void trace_critical_arclength_binary(const ModelType& model, double T0, const VecType &rhovec0, const std::string &filename) {
double t = 0.0, dt = 100;
- std::valarray<double> last_drhodt;
+ VecType last_drhodt;
VecType rhovec = rhovec0;
double T = T0;
@@ -299,17 +306,17 @@ void trace_critical_arclength_binary(const ModelType& model, double T0, const Ve
}
}
- auto drhodT = get_drhovec_dT_crit(model, T, rhovec);
+ auto drhodT = get_drhovec_dT_crit(model, T, rhovec).array().eval();
auto dTdt = 1.0 / norm(drhodT);
auto drhodt = drhodT * dTdt;
- auto eval = [](const auto& ex) { return std::valarray<bool>(ex); };
-
// Flip the sign if the tracing wants to go backwards, or if the first step would take you to negative concentrations
- if (iter == 0 && any(eval((rhovec + c * drhodt * dt) < 0))) {
+ Eigen::ArrayXd step = rhovec + c*drhodt*dt;
+ auto negativestepvals = (step < 0).eval();
+ if (iter == 0 && negativestepvals.all()) {
c *= -1;
}
- else if (iter > 0 && dot(std::valarray<double>(c * drhodt), last_drhodt) < 0) {
+ else if (iter > 0 && dot((c * drhodt).eval(), last_drhodt) < 0) {
c *= -1;
}
@@ -323,7 +330,7 @@ void trace_critical_arclength_binary(const ModelType& model, double T0, const Ve
try {
auto [Tnew, rhovecnew] = critical_polish_molefrac(model, T, rhovec, z0);
- T = Tnew; rhovec = VecType(&(rhovecnew[0]), rhovecnew.size());
+ T = Tnew; rhovec = rhovecnew;
}
catch (std::exception& e) {
std::cout << e.what() << std::endl;
diff --git a/include/teqp/derivs.hpp b/include/teqp/derivs.hpp
index ef9ed41..aad41a8 100644
--- a/include/teqp/derivs.hpp
+++ b/include/teqp/derivs.hpp
@@ -179,7 +179,7 @@ auto get_Bnvir(const Model& model, const TType &T, const ContainerType& molefrac
}
template <typename Model, typename TType, typename ContainerType>
-typename ContainerType::value_type get_B12vir(const Model& model, const TType T, const ContainerType& molefrac) {
+auto get_B12vir(const Model& model, const TType T, const ContainerType& molefrac) {
auto B2 = get_B2vir(model, T, molefrac); // Overall B2 for mixture
auto B20 = get_B2vir(model, T, std::valarray<double>({ 1,0 })); // Pure first component with index 0
@@ -199,15 +199,6 @@ typename ContainerType::value_type get_splus(const Model& model, const TType T,
return model.alphar(T, rhotot, molefrac) - get_Ar10(model, T, rhovec);
}
-/***
-* \brief Calculate Psir=ar*rho
-*/
-template <typename TType, typename ContainerType, typename Model>
-typename ContainerType::value_type get_Psir(const Model& model, const TType T, const ContainerType& rhovec) {
- auto rhotot_ = std::accumulate(std::begin(rhovec), std::end(rhovec), (decltype(rhovec[0]))0.0);
- return model.alphar(T, rhotot_, rhovec / rhotot_) * model.R * T * rhotot_;
-}
-
/***
* \brief Calculate the residual pressure from derivatives of alphar
*/
@@ -218,46 +209,97 @@ typename ContainerType::value_type get_pr(const Model& model, const TType T, con
return get_Ar01(model, T, rhotot_, rhovec / rhotot_) * rhotot_ * model.R * T;
}
+template<typename Model, typename Scalar=double, typename VectorType = Eigen::ArrayXd>
+struct IsochoricDerivatives{
+ /***
+ * \brief Calculate Psir=ar*rho
+ */
+ static auto get_Psir(const Model& model, const Scalar &T, const VectorType& rhovec) {
+ auto rhotot_ = std::accumulate(std::begin(rhovec), std::end(rhovec), (decltype(rhovec[0]))0.0);
+ return model.alphar(T, rhotot_, rhovec / rhotot_) * model.R * T * rhotot_;
+ }
-/***
-* \brief Calculate the Hessian of Psir = ar*rho w.r.t. the molar concentrations
-*
-* Requires the use of autodiff derivatives to calculate second partial derivatives
-*/
-template<typename Model, typename TType, typename RhoType>
-auto build_Psir_Hessian_autodiff(const Model& model, const TType &T, const RhoType& rho) {
- // Double derivatives in each component's concentration
- // N^N matrix (symmetric)
+ /***
+ * \brief Calculate the Hessian of Psir = ar*rho w.r.t. the molar concentrations
+ *
+ * Requires the use of autodiff derivatives to calculate second partial derivatives
+ */
+ static auto build_Psir_Hessian_autodiff(const Model& model, const Scalar& T, const VectorType& rho) {
+ // Double derivatives in each component's concentration
+ // N^N matrix (symmetric)
- dual2nd u; // the output scalar u = f(x), evaluated together with Hessian below
- VectorXdual2nd g;
- VectorXdual2nd rhovecc(rho.size()); for (auto i = 0; i < rho.size(); ++i) { rhovecc[i] = rho[i]; }
- auto hfunc = [&model, &T](const VectorXdual2nd& rho_) {
- auto rhotot_ = rho_.sum();
- auto molefrac = (rho_ / rhotot_).eval();
- return eval(model.alphar(T, rhotot_, molefrac) * model.R * T * rhotot_);
- };
- return autodiff::hessian(hfunc, wrt(rhovecc), at(rhovecc), u, g).eval(); // evaluate the function value u, its gradient, and its Hessian matrix H
-}
+ dual2nd u; // the output scalar u = f(x), evaluated together with Hessian below
+ VectorXdual2nd g;
+ VectorXdual2nd rhovecc(rho.size()); for (auto i = 0; i < rho.size(); ++i) { rhovecc[i] = rho[i]; }
+ auto hfunc = [&model, &T](const VectorXdual2nd& rho_) {
+ auto rhotot_ = rho_.sum();
+ auto molefrac = (rho_ / rhotot_).eval();
+ return eval(model.alphar(T, rhotot_, molefrac) * model.R * T * rhotot_);
+ };
+ return autodiff::hessian(hfunc, wrt(rhovecc), at(rhovecc), u, g).eval(); // evaluate the function value u, its gradient, and its Hessian matrix H
+ }
-/***
-* \brief Calculate the Hessian of Psir = ar*rho w.r.t. the molar concentrations
-*
-* Requires the use of multicomplex derivatives to calculate second partial derivatives
-*/
-template<typename Model, typename TType, typename RhoType>
-auto build_Psir_Hessian_mcx(const Model& model, const TType &T, const RhoType& rho) {
- // Double derivatives in each component's concentration
- // N^N matrix (symmetric)
- using namespace mcx;
+ /***
+ * \brief Calculate the Hessian of Psi = a*rho w.r.t. the molar concentrations
+ *
+ * Uses autodiff derivatives to calculate second partial derivatives
+ */
+ static auto build_Psi_Hessian_autodiff(const Model& model, const Scalar& T, const VectorType& rho) {
+ auto H = build_Psir_Hessian_autodiff(model, T, rho);
+ for (auto i = 0; i < 2; ++i) {
+ H(i, i) += model.R * T / rho[i];
+ }
+ return H;
+ }
+
+ /***
+ * \brief Calculate the Hessian of Psir = ar*rho w.r.t. the molar concentrations (residual contribution only)
+ *
+ * Requires the use of multicomplex derivatives to calculate second partial derivatives
+ */
+ static auto build_Psir_Hessian_mcx(const Model& model, const Scalar& T, const VectorType& rho) {
+ // Double derivatives in each component's concentration
+ // N^N matrix (symmetric)
+ using namespace mcx;
- // Lambda function for getting Psir with multicomplex concentrations
- using fcn_t = std::function< MultiComplex<double>(const std::valarray<MultiComplex<double>>&)>;
- fcn_t func = [&model, &T](const auto& rhovec) {
- return get_Psir(model, T, rhovec);
- };
- using mattype = Eigen::ArrayXXd;
- auto H = get_Hessian<mattype, fcn_t, std::valarray<double>, HessianMethods::Multiple>(func, rho);
- return H;
-}
\ No newline at end of file
+ // Lambda function for getting Psir with multicomplex concentrations
+ using fcn_t = std::function< MultiComplex<double>(const Eigen::ArrayX<MultiComplex<double>>&)>;
+ fcn_t func = [&model, &T](const auto& rhovec) {
+ auto rhotot_ = rhovec.sum();
+ auto molefrac = (rhovec / rhotot_).eval();
+ return model.alphar(T, rhotot_, molefrac) * model.R * T * rhotot_;
+ };
+ using mattype = Eigen::ArrayXXd;
+ auto H = get_Hessian<mattype, fcn_t, VectorType, HessianMethods::Multiple>(func, rho);
+ return H;
+ }
+
+ /***
+ * \brief Gradient of Psir = ar*rho w.r.t. the molar concentrations
+ *
+ * Uses autodiff to calculate second partial derivatives
+ */
+ static auto build_Psir_gradient_autodiff(const Model& model, const Scalar& T, const VectorType& rho) {
+ VectorXdual2nd rhovecc(rho.size()); for (auto i = 0; i < rho.size(); ++i) { rhovecc[i] = rho[i]; }
+ auto psirfunc = [&model, &T](const VectorXdual2nd& rho_) {
+ auto rhotot_ = rho_.sum();
+ auto molefrac = (rho_ / rhotot_).eval();
+ return eval(model.alphar(T, rhotot_, molefrac) * model.R * T * rhotot_);
+ };
+ auto val = autodiff::gradient(psirfunc, wrt(rhovecc), at(rhovecc)).eval(); // evaluate the gradient
+ return val;
+ }
+
+ /***
+ * \brief Calculate the chemical potential of each component
+ *
+ * Uses autodiff derivatives to calculate second partial derivatives
+ * See Eq. 9 of https://doi.org/10.1002/aic.16730
+ * \note: Some contributions to the ideal gas part are missing (reference state and cp0), but are not relevant to phase equilibria
+ */
+ static auto get_chempot_autodiff(const Model& model, const Scalar& T, const VectorType& rho) {
+ typename VectorType::value_type rhotot = rho.sum();
+ return (build_Psir_gradient_autodiff(model, T, rho).array() + model.R*T*(1.0 + log(rho / rhotot))).eval();
+ }
+};
\ No newline at end of file
diff --git a/src/multifluid.cpp b/src/multifluid.cpp
index 6b90279..5554926 100644
--- a/src/multifluid.cpp
+++ b/src/multifluid.cpp
@@ -38,7 +38,7 @@ void trace_critical_loci(const std::string &coolprop_root, const nlohmann::json
const auto &model = optmodel.value();
auto rhoc0 = 1.0 / model.redfunc.vc[i];
auto T0 = model.redfunc.Tc[i];
- std::valarray<double> rhovec(2); rhovec[i] = { rhoc0 }; rhovec[1L - i] = 0.0;
+ Eigen::ArrayXd rhovec(2); rhovec[i] = { rhoc0 }; rhovec[1L - i] = 0.0;
// Non-analytic terms make it impossible to initialize AT the pure components
if (pp[0] == "CarbonDioxide" || pp[1] == "CarbonDioxide") {
if (i == 0) {
@@ -70,18 +70,24 @@ int main(){
std::ifstream(coolprop_root + "/dev/mixtures/mixture_binary_pairs.json")
);
- // Critical curves
- {
- Timer t(1);
- trace_critical_loci(coolprop_root, BIPcollection);
- }
+ //// Critical curves
+ //{
+ // Timer t(1);
+ // trace_critical_loci(coolprop_root, BIPcollection);
+ //}
{
auto model = build_multifluid_model({ "methane", "ethane" }, coolprop_root, BIPcollection);
- std::valarray<double> rhovec = { 1.0, 2.0 };
+ Eigen::ArrayXd rhovec(2); rhovec << 1.0, 2.0 ;
double T = 300;
{
- const std::valarray<double> molefrac = { rhovec[0]/rhovec.sum(), rhovec[1]/rhovec.sum() };
+ const auto molefrac = (Eigen::ArrayXd(2) << rhovec[0]/rhovec.sum(), rhovec[1]/rhovec.sum()).finished();
+
+ auto B12 = get_B12vir(model, T, molefrac);
+
+ using id = IsochoricDerivatives<decltype(model)>;
+ auto mu = id::get_chempot_autodiff(model, T, rhovec);
+
const double rho = rhovec.sum();
volatile double T = 300.0;
constexpr int N = 10000;
@@ -93,7 +99,7 @@ int main(){
for (auto i = 0; i < N; ++i){
alphar = model.alphar(T, rho, molefrac);
}
- std::cout << alphar << std::endl;
+ std::cout << alphar << " function call" << std::endl;
}
{
Timer t(N);
@@ -121,25 +127,40 @@ int main(){
for (auto i = 0; i < N; ++i) {
alphar = get_Ar02(model, T, rho, molefrac);
}
- std::cout << alphar << std::endl;
+ std::cout << alphar << "; 2nd autodiff" << std::endl;
+ }
+ {
+ Timer t(N);
+ for (auto i = 0; i < N; ++i) {
+ auto o = 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 = 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 = get_Bnvir<5, ADBackends::autodiff>(model, T, molefrac)[5];
+ }
+ std::cout << alphar << "; 5 derivs" << std::endl;
}
}
+ const auto molefrac = (Eigen::ArrayXd(2) << 1.0 / 3.0, 2.0 / 3.0 ).finished();
+
auto alphar = model.alphar(T, rhovec);
auto Ar01 = get_Ar01(model, T, rhovec);
auto Ar10 = get_Ar10(model, T, rhovec);
+ auto Ar02 = get_Ar02(model, T, rhovec.sum(), molefrac);
auto splus = get_splus(model, T, rhovec);
-
- std::valarray<double> molefrac = { 1.0/3.0, 2.0/3.0 };
- auto B2 = get_B2vir(model, T, molefrac);
-
- std::valarray<double> dilrho = 0.00000000001*molefrac;
- auto B2other = get_Ar01(model, T, dilrho)/dilrho.sum();
-
- std::valarray<double> zerorho = 0.0*rhovec;
- auto Ar01dil = get_Ar01(model, T, zerorho);
- int ttt =0 ;
+ int ttt = 0;
}
return EXIT_SUCCESS;
}
diff --git a/src/tests/catch_tests.cxx b/src/tests/catch_tests.cxx
index 32377d5..0cdccf4 100644
--- a/src/tests/catch_tests.cxx
+++ b/src/tests/catch_tests.cxx
@@ -84,12 +84,14 @@ TEST_CASE("Check virial coefficients for vdW", "[virial]")
TEST_CASE("Check Hessian of Psir", "[virial]")
{
+
double T = 298.15;
double rho = 3.0;
- const std::valarray<double> rhovec = { rho / 2, rho / 2 };
+ const Eigen::Array2d rhovec = { rho / 2, rho / 2 };
auto get_worst_error = [&T, &rhovec](const auto &model){
- auto H1 = build_Psir_Hessian_autodiff(model, T, rhovec);
- auto H2 = build_Psir_Hessian_mcx(model, T, rhovec);
+ using id = IsochoricDerivatives <decltype(model)>;
+ auto H1 = id::build_Psir_Hessian_autodiff(model, T, rhovec);
+ auto H2 = id::build_Psir_Hessian_mcx(model, T, rhovec);
auto err = (H1.array() - H2).abs().maxCoeff();
CAPTURE(err);
CHECK(err < 1e-15);
@@ -139,7 +141,7 @@ TEST_CASE("Trace critical locus for vdW", "[vdW][crit]")
auto Zc = 3.0/8.0;
auto rhoc0 = pc_Pa[0] / (vdW.R * Tc_K[0]) / Zc;
double T0 = Tc_K[0];
- std::valarray<double> rhovec0 = { rhoc0, 0.0 };
+ Eigen::ArrayXd rhovec0(2); rhovec0 << rhoc0, 0.0 ;
auto tic0 = std::chrono::steady_clock::now();
std::string filename = "";
@@ -152,4 +154,27 @@ TEST_CASE("TEST B12", "") {
const double T = 298.15;
const std::valarray<double> molefrac = { 1/3, 2/3 };
auto B12 = get_B12vir(model, T, molefrac);
+}
+
+TEST_CASE("Test psir gradient", "") {
+ const auto model = build_vdW();
+ const double T = 298.15;
+
+ using id = IsochoricDerivatives<decltype(model)>;
+ const Eigen::Array2d rhovec = { 1, 2 };
+ const Eigen::Array2d molefrac = { 1 / 3, 2 / 3 };
+ auto psirfunc2 = [&model](const auto& T, const auto& rho_) {
+ auto rhotot_ = rho_.sum();
+ auto molefrac = (rho_ / rhotot_);
+ return model.alphar(T, rhotot_, molefrac) * model.R * T * rhotot_;
+ };
+ auto chk0 = derivrhoi(psirfunc2, T, rhovec, 0);
+ auto chk1 = derivrhoi(psirfunc2, T, rhovec, 1);
+ auto grad = id::build_Psir_gradient_autodiff(model, T, rhovec);
+ auto err0 = std::abs((chk0 - grad[0])/chk0);
+ auto err1 = std::abs((chk1 - grad[1])/chk1);
+ CAPTURE(err0);
+ CAPTURE(err1);
+ CHECK(err0 < 1e-12);
+ CHECK(err1 < 1e-12);
}
\ No newline at end of file
--
GitLab