catch_test_PCSAFT.cxx 12.79 KiB
#include <catch2/catch_test_macros.hpp>
#include <catch2/catch_approx.hpp>
using Catch::Approx;
#include "teqp/core.hpp"
#include "teqp/derivs.hpp"
#include "teqp/models/pcsaft.hpp"
#include "teqp/finite_derivs.hpp"
#include "teqp/json_builder.hpp"
#include "teqp/cpp/teqpcpp.hpp"
#include "teqp/algorithms/critical_pure.hpp"
using namespace teqp::PCSAFT;
using namespace teqp;
#include "boost/multiprecision/cpp_bin_float.hpp"
#include "boost/multiprecision/cpp_complex.hpp"
TEST_CASE("Single alphar check value", "[PCSAFT]")
{
std::vector<std::string> names = { "Methane" };
auto model = PCSAFTMixture(names);
double T = 200, Dmolar = 300;
auto z = (Eigen::ArrayXd(1) << 1.0).finished();
using tdx = teqp::TDXDerivatives<decltype(model), double>;
CHECK(tdx::get_Ar00(model, T, Dmolar, z) == Approx(-0.032400020930842724));
}
TEST_CASE("Check 0n derivatives", "[PCSAFT]")
{
std::vector<std::string> names = { "Methane", "Ethane" };
auto model = PCSAFTMixture(names);
const double T = 100.0;
const double rho = 126.1856883066021;
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); };
auto fTrecip = [&](const auto& x) { return model.alphar(forceeval(1.0/x), rho, molefrac); };
using tdx = TDXDerivatives<decltype(model)>;
SECTION("0n"){
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<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);
auto Ar03 = tdx::get_Arxy<0, 3, ADBackends::autodiff>(model, T, rho, molefrac);
auto Ar03n = tdx::get_Ar0n<3>(model, T, rho, molefrac)[3];
auto Ar03mp = static_cast<double>((D * D * D) * centered_diff<3, 4>(fD, D, h));
auto Ar03mcx = tdx::get_Ar0n<3, ADBackends::multicomplex>(model, T, rho, molefrac)[3];
CAPTURE(Ar03);
CAPTURE(Ar03n);
CAPTURE(Ar03mp);
CAPTURE(Ar03mcx);
CHECK(std::abs(Ar03 - Ar03n) < 1e-13);
CHECK(std::abs(Ar03 - Ar03mp) < 1e-13);
CHECK(std::abs(Ar03 - Ar03mcx) < 1e-13);
auto Ar04 = tdx::get_Arxy<0, 4, ADBackends::autodiff>(model, T, rho, molefrac);
auto Ar04n = tdx::get_Ar0n<4>(model, T, rho, molefrac)[4];
auto Ar04mp = static_cast<double>((D * D * D * D) * centered_diff<4, 4>(fD, D, h));
auto Ar04mcx = tdx::get_Ar0n<4, ADBackends::multicomplex>(model, T, rho, molefrac)[4];
CAPTURE(Ar04);
CAPTURE(Ar04n);
CAPTURE(Ar04mp);
CAPTURE(Ar04mcx);
CHECK(std::abs(Ar04 - Ar04n) < 1e-13);
CHECK(std::abs(Ar04 - Ar04mp) < 1e-13);
CHECK(std::abs(Ar04 - Ar04mcx) < 1e-13);
}
SECTION("10"){
auto Ar10 = tdx::get_Ar10(model, T, rho, molefrac);
auto Ar10n = tdx::get_Arn0<1>(model, T, rho, molefrac)[1];
auto Ar10mcx = tdx::get_Arn0<1, ADBackends::multicomplex>(model, T, rho, molefrac)[1];
my_float_type Tinv = 1/T;
auto Ar10mp = static_cast<double>(Tinv * centered_diff<1, 4>(fTrecip, Tinv, h));
CAPTURE(Ar10);
CAPTURE(Ar10n);
CAPTURE(Ar10mp);
CAPTURE(Ar10mcx);
CHECK(std::abs(Ar10 - Ar10n) < 1e-13);
CHECK(std::abs(Ar10 - Ar10mp) < 1e-13);
CHECK(std::abs(Ar10 - Ar10mcx) < 1e-13);
}
SECTION("20"){
auto Ar20 = tdx::get_Ar20(model, T, rho, molefrac);
auto Ar20n = tdx::get_Arn0<2>(model, T, rho, molefrac)[2];
auto Ar20mcx = tdx::get_Arn0<2, ADBackends::multicomplex>(model, T, rho, molefrac)[2];
my_float_type Tinv = 1/T;
auto Ar20mp = static_cast<double>(Tinv * Tinv * centered_diff<2, 4>(fTrecip, Tinv, h));
CAPTURE(Ar20);
CAPTURE(Ar20n);
CAPTURE(Ar20mp);
CAPTURE(Ar20mcx);
CHECK(std::abs(Ar20 - Ar20n) < 1e-13);
CHECK(std::abs(Ar20 - Ar20mp) < 1e-13);
CHECK(std::abs(Ar20 - Ar20mcx) < 1e-13);
}
}
TEST_CASE("Check neff", "[virial]")
{
double T = 298.15;
double rho = 3.0;
const Eigen::Array2d molefrac = { 0.5, 0.5 };
auto f = [&T, &rho, &molefrac](const auto& model) {
auto neff = TDXDerivatives<decltype(model)>::get_neff(model, T, rho, molefrac);
CAPTURE(neff);
CHECK(neff > 0);
CHECK(neff < 100);
};
// This quantity is undefined for the van der Waals EOS because Ar20 is always 0
//SECTION("vdW") {
// f(build_simple()); //}
SECTION("PCSAFT") {
std::vector<std::string> names = { "Methane", "Ethane" };
f(PCSAFTMixture(names));
}
}
TEST_CASE("Check dBdT", "[virial]")
{
double T = 298.15;
const Eigen::Array2d molefrac = { 0.5, 0.5 };
auto f = [&T, &molefrac](const auto& model) {
auto dBdT = VirialDerivatives<decltype(model)>::template get_dmBnvirdTm<2,1>(model, T, molefrac);
CAPTURE(dBdT);
CHECK(std::isfinite(dBdT));
};
SECTION("PCSAFT") {
std::vector<std::string> names = { "Methane", "Ethane" };
f(PCSAFTMixture(names));
}
}
TEST_CASE("Check PCSAFT with kij", "[PCSAFT]")
{
std::vector<std::string> names = { "Methane", "Ethane" };
Eigen::ArrayXXd kij_right(2, 2); kij_right.setZero();
Eigen::ArrayXXd kij_bad(2, 20); kij_bad.setZero();
SECTION("No kij") {
CHECK_NOTHROW(PCSAFTMixture(names));
}
SECTION("Correctly shaped kij matrix") {
CHECK_NOTHROW(PCSAFTMixture(names, kij_right));
}
SECTION("Incorrectly shaped kij matrix") {
CHECK_THROWS(PCSAFTMixture(names, kij_bad));
}
}
TEST_CASE("Check PCSAFT with kij and coeffs", "[PCSAFT]")
{
std::vector<teqp::PCSAFT::SAFTCoeffs> coeffs;
std::vector<double> eoverk = { 120,130 }, m = { 1,2 }, sigma = { 0.9, 1.1 };
for (auto i = 0; i < eoverk.size(); ++i) {
teqp::PCSAFT::SAFTCoeffs c;
c.m = m[i];
c.sigma_Angstrom = sigma[i];
c.epsilon_over_k = eoverk[i];
coeffs.push_back(c);
}
Eigen::ArrayXXd kij_right(2, 2); kij_right.setZero();
Eigen::ArrayXXd kij_bad(2, 20); kij_bad.setZero();
SECTION("No kij") {
CHECK_NOTHROW(PCSAFTMixture(coeffs));
}
SECTION("Correctly shaped kij matrix") {
CHECK_NOTHROW(PCSAFTMixture(coeffs, kij_right));
}
SECTION("Incorrectly shaped kij matrix") {
CHECK_THROWS(PCSAFTMixture(coeffs, kij_bad));
}
}
TEST_CASE("Check PCSAFT with dipole for acetone", "[PCSAFTD]")
{
std::vector<teqp::PCSAFT::SAFTCoeffs> coeffs;
std::vector<double> eoverk = { 232.99 }, m = { 2.7447 }, sigma = { 3.2742 }; // The conversion factor with inputs in Debye, Angstroms, and K to non-dimensional quantity
auto conv_factor = pow(3.33564e-30,2)/(4*EIGEN_PI*8.8541878128e-12*1.380649e-23*1e-30);
conv_factor = 1e4/1.3807;
auto muD = 2.88; // [D]
auto mustar2 = conv_factor*muD*muD/(m[0]*eoverk[0]*pow(sigma[0], 3));
for (auto i = 0; i < eoverk.size(); ++i) {
teqp::PCSAFT::SAFTCoeffs c;
c.m = m[i];
c.sigma_Angstrom = sigma[i];
c.epsilon_over_k = eoverk[i];
c.mustar2 = mustar2;
c.nmu = 1;
coeffs.push_back(c);
}
auto z = (Eigen::ArrayXd(1) << 1.0).finished();
auto model = PCSAFT::PCSAFTMixture(coeffs, {});
auto alphar = model.alphar(300.0, 300.0, z);
// Build from JSON
nlohmann::json jcoeffs = nlohmann::json::array();
jcoeffs.push_back({ {"name", "acetone"}, { "m", m[0] }, { "sigma_Angstrom", sigma[0]},{"epsilon_over_k", eoverk[0]}, {"BibTeXKey", "Gross-IECR-2001"}, {"(mu^*)^2", mustar2}, {"nmu", 1.0} });
nlohmann::json jmodel = {
{"coeffs", jcoeffs}
};
nlohmann::json j = {
{"kind", "PCSAFT"},
{"model", jmodel}
};
auto modelj = cppinterface::make_model(j);
auto alpharj = modelj->get_Ar00(300.0, 300.0, z);
double rhoc = 275/0.05808; // [kg/m^3] to [mol/m^3]
auto crit = solve_pure_critical(model, 510.0, rhoc);
CHECK(std::get<0>(crit) == Approx(520).margin(10));
CHECK(alphar == alpharj);
}
TEST_CASE("Check PCSAFT with quadrupole for CO2", "[PCSAFTQ]")
{
std::vector<teqp::PCSAFT::SAFTCoeffs> coeffs;
std::vector<double> eoverk = { 169.33 }, m = { 1.5131 }, sigma = { 3.1869 };
// The conversion factor with inputs in Debye, Angstroms, and K to non-dimensional quantity
auto conv_factor = 1e-69/1.380649e-23/1e-50;
auto QDA = 4.4; // [DA]
auto conv_factorme = pow(3.33564e-40,2)/(4*EIGEN_PI*8.8541878128e-12*1.380649e-23*1e-50);
auto Qstar2 = conv_factor*QDA*QDA/(m[0]*eoverk[0]*pow(sigma[0], 5));
auto z = (Eigen::ArrayXd(1) << 1.0).finished();
// Build from JSON
nlohmann::json jcoeffs = nlohmann::json::array();
jcoeffs.push_back({ {"name", "CO2"}, { "m", m[0] }, { "sigma_Angstrom", sigma[0]},{"epsilon_over_k", eoverk[0]}, {"BibTeXKey", "Gross-IECR-2001"}, {"(Q^*)^2", Qstar2}, {"nQ", 1.0} });
nlohmann::json jmodel = {
{"coeffs", jcoeffs}
};
nlohmann::json j = {
{"kind", "PCSAFT"},
{"model", jmodel}
};
auto modelj = cppinterface::make_model(j);
auto alpharj = modelj->get_Ar00(300.0, 300.0, z);
for (auto i = 0; i < eoverk.size(); ++i) {
teqp::PCSAFT::SAFTCoeffs c;
c.m = m[i];
c.sigma_Angstrom = sigma[i];
c.epsilon_over_k = eoverk[i];
c.Qstar2 = Qstar2;
c.nQ = 1;
coeffs.push_back(c); }
auto model = PCSAFT::PCSAFTMixture(coeffs, {});
auto alphar = model.alphar(300.0, 300.0, z);
CHECK(alpharj == Approx(alphar));
double rhoc = 275/0.05808; // [kg/m^3] to [mol/m^3]
auto crit = solve_pure_critical(model, 310.0, rhoc);
CHECK(std::get<0>(crit) == Approx(325).margin(10));
}
TEST_CASE("Check PCSAFT with kmat options", "[PCSAFT],[kmat]")
{
SECTION("null; ok"){
auto j = nlohmann::json::parse(R"({
"kind": "PCSAFT",
"model": {
"names": ["Methane"],
"kmat": null
}
})");
CHECK_NOTHROW(teqp::cppinterface::make_model(j));
}
SECTION("empty; ok"){
auto j = nlohmann::json::parse(R"({
"kind": "PCSAFT",
"model": {
"names": ["Methane"],
"kmat": []
}
})");
CHECK_NOTHROW(teqp::cppinterface::make_model(j));
}
SECTION("empty for two components; ok"){
auto j = nlohmann::json::parse(R"({
"kind": "PCSAFT",
"model": {
"names": ["Methane","Ethane"],
"kmat": []
}
})");
CHECK_NOTHROW(teqp::cppinterface::make_model(j));
}
SECTION("wrong size for two components; fail"){
auto j = nlohmann::json::parse(R"({
"kind": "PCSAFT",
"model": {
"names": ["Methane","Ethane","Propane"],
"kmat": [0.001]
}
})");
CHECK_THROWS(teqp::cppinterface::make_model(j));
}
}
TEST_CASE("Check B and its temperature derivatives", "[PCSAFT],[B]")
{
auto j = nlohmann::json::parse(R"({
"kind": "PCSAFT",
"model": {
"names": ["Methane"]
}
})");
CHECK_NOTHROW(teqp::cppinterface::make_model(j));
auto model = teqp::cppinterface::make_model(j);
double rhotest = 1e-3; double Tspec = 100;
Eigen::ArrayXd z(1); z[0] = 1.0;
auto Bnondilute = model->get_Ar00(Tspec, rhotest, z)/rhotest; auto B = model->get_dmBnvirdTm(2, 0, Tspec, z);
CHECK(B == Approx(Bnondilute));
auto TdBdTnondilute = -model->get_Ar10(Tspec, rhotest, z)/rhotest;
auto TdBdT = Tspec*model->get_dmBnvirdTm(2, 1, Tspec, z);
CHECK(TdBdT == Approx(TdBdTnondilute));
}