From c408d38815e9274fd0d1ab6801ce51f72d674b27 Mon Sep 17 00:00:00 2001 From: Ian Bell <ian.bell@nist.gov> Date: Sun, 31 Jul 2022 15:42:37 -0400 Subject: [PATCH] And this time with a test --- src/tests/catch_test_cubics.cxx | 56 +++++++++++++++++++++++++++++++++ 1 file changed, 56 insertions(+) diff --git a/src/tests/catch_test_cubics.cxx b/src/tests/catch_test_cubics.cxx index ef98d14..838d253 100644 --- a/src/tests/catch_test_cubics.cxx +++ b/src/tests/catch_test_cubics.cxx @@ -185,4 +185,60 @@ TEST_CASE("Check manual integration of subcritical VLE isotherm for binary mixtu double pfinal = J.back().at("pL / Pa").back(); CHECK(std::abs(pfinal / pfinal_goal-1) < 1e-5); } +} + +TEST_CASE("Check infinite dilution of subcritical VLE derivatives", "[cubic][isochoric][infdil]") +{ + // Values taken from http://dx.doi.org/10.6028/jres.121.011 + std::valarray<double> Tc_K = { 190.564, 154.581 }, + pc_Pa = { 4599200, 5042800 }, + acentric = { 0.011, 0.022 }; + auto model = canonical_PR(Tc_K, pc_Pa, acentric); + const auto N = Tc_K.size(); + using state_type = std::valarray<double>; + REQUIRE(N == 2); + auto get_start = [&](double T, auto i) { + std::valarray<double> Tc_(Tc_K[i], 1), pc_(pc_Pa[i], 1), acentric_(acentric[i], 1); + auto PR = canonical_PR(Tc_, pc_, acentric_); + auto [rhoL, rhoV] = PR.superanc_rhoLV(T); + state_type o(N * 2); + o[i] = rhoL; + o[i + N] = rhoV; + return o; + }; + int i = 1; + double T = 120; + auto rhostart_dil = get_start(T, i); + auto rhostart_notdil = rhostart_dil; + rhostart_notdil[1-i] += 1e-6; + rhostart_notdil[1-i+N] += 1e-6; + + SECTION("Along isotherm") { + // Derivative function with respect to p + auto xprime = [&](const state_type& X) { + REQUIRE(X.size() % 2 == 0); + auto N = X.size() / 2; + // Memory maps into the state vector for inputs and their derivatives + auto rhovecL = Eigen::Map<const Eigen::ArrayXd>(&(X[0]), N); + auto rhovecV = Eigen::Map<const Eigen::ArrayXd>(&(X[0]) + N, N); + return get_drhovecdp_Tsat(model, T, rhovecL, rhovecV); + }; + auto dernotdil = xprime(rhostart_notdil); + auto derdil = xprime(rhostart_dil); + CHECK(true); + } + SECTION("Along isobar") { + // Derivative function with respect to T + auto xprime = [&](const state_type& X) { + REQUIRE(X.size() % 2 == 0); + auto N = X.size() / 2; + // Memory maps into the state vector for inputs and their derivatives + auto rhovecL = Eigen::Map<const Eigen::ArrayXd>(&(X[0]), N); + auto rhovecV = Eigen::Map<const Eigen::ArrayXd>(&(X[0]) + N, N); + return get_drhovecdT_psat(model, T, rhovecL.eval(), rhovecV.eval()); + }; + auto dernotdil = xprime(rhostart_notdil); + auto derdil = xprime(rhostart_dil); + CHECK(true); + } } \ No newline at end of file -- GitLab