Skip to content

Commit

Permalink
Move PC-SAFT catch tests into own file
Browse files Browse the repository at this point in the history 
Should speed up compilation a smidge
  • Loading branch information
@ianhbell
ianhbell committed Aug 10, 2022
1 parent d406028 commit aa09674
Show file tree
Hide file tree
Showing 2 changed files with 92 additions and 89 deletions.
92 changes: 92 additions & 0 deletions src/tests/catch_test_PCSAFT.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,11 @@
using Catch::Approx;

#include "teqp/core.hpp"
#include "teqp/derivs.hpp"
#include "teqp/models/pcsaft.hpp"
#include "teqp/finite_derivs.hpp"
using namespace teqp::PCSAFT;
using namespace teqp;

TEST_CASE("Single alphar check value", "[PCSAFT]")
{
Expand All @@ -17,6 +20,95 @@ TEST_CASE("Single alphar check value", "[PCSAFT]")
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); };

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<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);
}

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 PCSAFT with kij", "[PCSAFT]")
{
std::vector<std::string> names = { "Methane", "Ethane" };
Expand Down
89 changes: 0 additions & 89 deletions src/tests/catch_tests.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,6 @@
using Catch::Approx;

#include "teqp/core.hpp"
#include "teqp/models/pcsaft.hpp"
using namespace teqp::PCSAFT;
#include "teqp/models/cubicsuperancillary.hpp"
#include "teqp/models/CPA.hpp"
#include "teqp/models/vdW.hpp"
Expand Down Expand Up @@ -100,26 +98,6 @@ TEST_CASE("Check virial coefficients for vdW", "[virial]")
}
}

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 Hessian of Psir", "[virial]")
{
Expand Down Expand Up @@ -177,73 +155,6 @@ TEST_CASE("Check p four ways for vdW", "[virial][p]")
CHECK(std::abs(p_ar0n - pexact) / pexact < 1e-8);
}

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); };

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<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);
}

TEST_CASE("Test all vdW derivatives good to numerical precision", "[vdW]")
{
using my_float_type = boost::multiprecision::number<boost::multiprecision::cpp_bin_float<100U>>;
Expand Down

0 comments on commit aa09674

Please sign in to comment.