From 06f046c533c186469a261c00b495fff03026b8a2 Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.bell@nist.gov>
Date: Fri, 27 May 2022 10:03:29 -0400
Subject: [PATCH] Extend the PC-SAFT derivatives
---
src/tests/catch_tests.cxx | 28 +++++++++++++++++++++++++---
1 file changed, 25 insertions(+), 3 deletions(-)
diff --git a/src/tests/catch_tests.cxx b/src/tests/catch_tests.cxx
index 97aee71..7d92a0f 100644
--- a/src/tests/catch_tests.cxx
+++ b/src/tests/catch_tests.cxx
@@ -177,7 +177,7 @@ TEST_CASE("Check p four ways for vdW", "[virial][p]")
CHECK(std::abs(p_ar0n - pexact) / pexact < 1e-8);
}
-TEST_CASE("Check 0n derivatives", "[virial][p]")
+TEST_CASE("Check 0n derivatives", "[PCSAFT]")
{
std::vector<std::string> names = { "Methane", "Ethane" };
auto model = PCSAFTMixture(names);
@@ -187,15 +187,37 @@ TEST_CASE("Check 0n derivatives", "[virial][p]")
const auto rhovec = (Eigen::ArrayXd(2) << rho, 0).finished();
const auto molefrac = rhovec / rhovec.sum();
+ using my_float_type = boost::multiprecision::number<boost::multiprecision::cpp_bin_float<100U>>;
+ my_float_type D = rho, h = pow(my_float_type(10.0), -10);
+ auto fD = [&](const auto& x) { return model.alphar(T, x, molefrac); };
+
using tdx = TDXDerivatives<decltype(model)>;
auto Ar02 = tdx::get_Ar02(model, T, rho, molefrac);
auto Ar02n = tdx::get_Ar0n<2>(model, T, rho, molefrac)[2];
+ auto Ar02mp = static_cast<double>((D*D)*centered_diff<2, 4>(fD, D, h));
+ auto Ar02mcx = tdx::get_Ar0n<2,ADBackends::multicomplex>(model, T, rho, molefrac)[2];
+ CAPTURE(Ar02);
+ CAPTURE(Ar02n);
+ CAPTURE(Ar02mp);
+ CAPTURE(Ar02mcx);
CHECK(std::abs(Ar02 - Ar02n) < 1e-13);
+ CHECK(std::abs(Ar02 - Ar02mp) < 1e-13);
+ CHECK(std::abs(Ar02 - Ar02mcx) < 1e-13);
auto Ar01 = tdx::get_Ar01(model, T, rho, molefrac);
- auto Ar01n = tdx::get_Ar0n<4>(model, T, rho, molefrac)[1];
+ auto Ar01n = tdx::get_Ar0n<1>(model, T, rho, molefrac)[1];
+ auto Ar01mcx = tdx::get_Ar0n<1,ADBackends::multicomplex>(model, T, rho, molefrac)[1];
+ auto Ar01csd = tdx::get_Ar01<ADBackends::complex_step>(model, T, rho, molefrac);
+ auto Ar01mp = static_cast<double>(D*centered_diff<1, 4>(fD, D, h));
+ CAPTURE(Ar01);
+ CAPTURE(Ar01n);
+ CAPTURE(Ar01mp);
+ CAPTURE(Ar01mcx);
+ CAPTURE(Ar01csd);
CHECK(std::abs(Ar01 - Ar01n) < 1e-13);
-
+ CHECK(std::abs(Ar01 - Ar01mp) < 1e-13);
+ CHECK(std::abs(Ar01 - Ar01mcx) < 1e-13);
+ CHECK(std::abs(Ar01 - Ar01csd) < 1e-13);
}
TEST_CASE("Test all vdW derivatives good to numerical precision", "[vdW]")
--
GitLab