sat_Z_accuracy.cpp 11.02 KiB
/***
* A script for testing the loss in precision of autodiff differentiation and comparing to the lost
* precision in REFPROP
*/
// Only this file gets the implementation
#define REFPROP_IMPLEMENTATION
#define REFPROP_FUNCTION_MODIFIER
#include "REFPROP_lib.h"
#undef REFPROP_FUNCTION_MODIFIER
#undef REFPROP_IMPLEMENTATION
#include <iostream>
#include <valarray>
#include "teqp/models/multifluid.hpp"
#include "teqp/derivs.hpp"
// Imports from boost
#include <boost/multiprecision/cpp_bin_float.hpp>
using namespace boost::multiprecision;
#include "teqp/finite_derivs.hpp"
/// A standalone implementation to be more in control of type promotion.
/// In the end this standalone implementation gives the same answer
/// This is for propane
template<typename T, typename Tau, typename Delta>
auto alphar_Lemmon2009(Tau tau, Delta delta)
{
const static std::valarray<T> d = { 4.0,1.0,1.0,2.0,2.0,1.0,3.0,6.0,6.0,2.0,3.0,1.0,1.0,1.0,2.0,2.0,4.0,1.0 },
n = { 0.042910051,1.7313671,-2.4516524,0.34157466,-0.46047898,-0.66847295,0.20889705,0.19421381,-0.22917851,-0.60405866,0.066680654,0.017534618,0.33874242,0.22228777,-0.23219062,-0.09220694,-0.47575718,-0.017486824 },
t = { 1,0.33,0.8,0.43,0.9,2.46,2.09,0.88,1.09,3.25,4.62,0.76,2.5,2.75,3.05,2.55,8.4,6.75 },
ld = { 0,0,0,0,0,1,1,1,1,2,2,0,0,0,0,0,0,0 },
cd = { 0,0,0,0,0,1,1,1,1,1,1,0,0,0,0,0,0,0 },
lt = { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
ct = { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
beta = { 0,0,0,0,0,0,0,0,0,0,0,2.33,3.47,3.15,3.19,0.92,18.8,547.8 },
epsilon = { 0,0,0,0,0,0,0,0,0,0,0,1.283,0.6936,0.788,0.473,0.8577,0.271,0.948 },
eta = { 0,0,0,0,0,0,0,0,0,0,0,0.963,1.977,1.917,2.307,2.546,3.28,14.6 },
gamma = { 0,0,0,0,0,0,0,0,0,0,0,0.684,0.829,1.419,0.817,1.5,1.426,1.093 };
std::common_type_t<Tau, Delta> result = 0.0;
for (auto i = 0; i < n.size(); ++i) {
result += n[i] * pow(tau, t[i]) * pow(delta, d[i]) * exp(-cd[i] * pow(delta, ld[i]) - eta[i] * pow(delta - epsilon[i], 2) - beta[i] * pow(tau - gamma[i], 2));
}
return result;
}
int REFPROP_setup(const std::string &RPname) {
// You may need to change this path to suit your installation
// Note: forward-slashes are recommended.
std::string path = std::getenv("RPPREFIX");
std::string DLL_name = "";
// Load the shared library and set up the fluid
std::string err;
bool loaded_REFPROP = load_REFPROP(err, path, DLL_name);
printf("Loaded refprop: %s @ address %zu\n", loaded_REFPROP ? "true" : "false", REFPROP_address());
if (!loaded_REFPROP) { throw std::invalid_argument("Bad load of REFPROP"); }
char hpath[256]; strcpy(hpath, (path + std::string(254-path.size(),'\0')).c_str());
SETPATHdll(hpath, 255);
int ierr = 0, nc = 1;
char herr[256], hfld[10000] = " ", hhmx[256] = "HMX.BNC", href[4] = "DEF";
strcpy(hfld, RPname.c_str());
SETUPdll(nc, hfld, hhmx, href, ierr, herr, 10000, 255, 3, 255);
if (ierr != 0) { throw std::invalid_argument("Bad setup of REFPROP: "+std::string(herr)); }
return 0;
}
struct REFPROP_sat_output {
double T, rhoLmol_L, rhoVmol_L, p_kPa, rho_mol_L; char herr[256];
int ierr;
std::valarray<double> mole_fractions{ 0.0, 20 }, mole_fractions_liq{ 0.0, 20 }, mole_fractions_vap{ 0.0, 20 };
};
// Do a saturation call in REFPROP to generate the liquid and vapor densities for a given temperature
auto REFPROP_sat(double T) {
REFPROP_sat_output o;
o.T = T;
int iFlsh = 0;
SATTdll(T, &(o.mole_fractions[0]), iFlsh,
o.p_kPa, o.rhoLmol_L, o.rhoVmol_L, &(o.mole_fractions_liq[0]), &(o.mole_fractions_vap[0]),
o.ierr, o.herr, errormessagelength);
return o;
}
struct calc_output {
double Zexact, Zteqp, Ar01exact, Ar01teqp, Ar02exact, Ar02teqp, Ar03exact, Ar03teqp;
};
template<typename Model, typename VECTOR>
auto with_teqp_and_boost(const Model &model, double T, double rho, const VECTOR &z, bool is_propane){
// Pressure for each phase via teqp in double precision w/ autodiff
using tdx = TDXDerivatives<decltype(model), double, VECTOR>;
double Zteqp = 1.0 + tdx::get_Ar01(model, T, rho, z);
double Ar01teqp = tdx::get_Ar01(model, T, rho, z);
// Calculation with ridiculous number of digits of precision (the approximation of ground truth)
using my_float = boost::multiprecision::number<boost::multiprecision::cpp_bin_float<200>>;
my_float Tc = model.redfunc.Tc[0];
my_float rhoc = 1.0/static_cast<my_float>(model.redfunc.vc[0]);
auto delta = static_cast<my_float>(rho) / rhoc;
auto tau = Tc / static_cast<my_float>(T);
my_float ddelta = 1e-30 * delta;
my_float deltaplus = delta + ddelta, deltaminus = delta - ddelta;
using coef_type = my_float; // What numerical type to use to initialize the coefficients (in the end it doesn't matter since they all get upcasted to my_float)
// Check that the function values are exactly the same
auto ar1 = model.corr.get_EOS(0).alphar(tau, delta);
if (is_propane) {
// As the standalone (if we are using propane)
auto ar2 = alphar_Lemmon2009<my_float>(tau, delta);
auto dar2 = static_cast<double>((ar2 - ar1) / ar1);
if (std::abs(dar2) > 1e-100) { // yes, we have ridiculously accurate values
throw std::invalid_argument("Function values are not exactly the same");
}
}
else {
// Or from REFPROP otherwise
int itau = 0, idelta = 0;
double tau_ = static_cast<double>(tau), delta_ = static_cast<double>(delta);
std::valarray<double> z(20); z = 1;
double ar2 = -1; PHIXdll(itau, idelta, tau_, delta_, &(z[0]), ar2);
double dar2 = static_cast<double>((ar2 - ar1) / ar1);
if (std::abs(dar2) > 5e-14) { // basically double precision..
std::cout << dar2 << std::endl;
throw std::invalid_argument("Function values are not exactly the same; error (%): "+std::to_string(dar2));
}
}
// And now the derivative value in two subtly different approaches, also check that 2nd-order-truncation and 4th-order-truncation are the same
auto derL2_2nd = (alphar_Lemmon2009<coef_type>(tau, deltaplus) - alphar_Lemmon2009<coef_type>(tau, deltaminus)) / (2.0 * ddelta) * delta;
auto derL2_4th = (
1.0*alphar_Lemmon2009<coef_type>(tau, delta - 2.0*ddelta)/12.0
-2.0*alphar_Lemmon2009<coef_type>(tau, delta - ddelta)/3.0
+2.0*alphar_Lemmon2009<coef_type>(tau, delta + ddelta)/3.0
-1.0*alphar_Lemmon2009<coef_type>(tau, delta + 2.0*ddelta)/12.0
) / ddelta * delta;
auto derL3 = (model.corr.get_EOS(0).alphar(tau, deltaplus) - model.corr.get_EOS(0).alphar(tau, deltaminus)) / (2.0 * ddelta) * delta;
auto Zexact = derL3 + 1.0;
if (is_propane) {
auto d3 = static_cast<double>((derL2_2nd - derL3) / derL2_2nd);
auto d34th = static_cast<double>((derL2_4th - derL2_2nd) / derL2_2nd);
if (std::abs(d3) > 1e-100) { // yes, we have ridiculously accurate values
throw std::invalid_argument("Derivatives are not exactly the same in teqp and in standalone implementation");
}
}
calc_output o;
o.Zexact = static_cast<double>(Zexact);
o.Zteqp = Zteqp;
o.Ar01exact = static_cast<double>(derL3);
o.Ar01teqp = Ar01teqp;
// Now do the third-order derivative of alphar, as a further test
// Define a generic lambda function taking rho
auto ff = [&](const auto& rho){ return model.alphar(T, rho, z); };
my_float drho = 1e-30*rho;
o.Ar02exact = static_cast<double>(centered_diff<2,6>(ff,static_cast<my_float>(rho),drho)*pow(rho, 2));
o.Ar02teqp = tdx::template get_Ar0n<2>(model, T, rho, z)[2];
o.Ar03exact = static_cast<double>(centered_diff<3,6>(ff,static_cast<my_float>(rho),drho)*pow(rho, 3));
o.Ar03teqp = tdx::template get_Ar0n<3>(model, T, rho, z)[3];
return o;
}
int do_one(const std::string &RPname, const std::string &teqpname)
{
REFPROP_setup(RPname);
auto model = build_multifluid_model({ teqpname }, "../mycp", "../mycp/dev/mixtures/mixture_binary_pairs.json");
Eigen::ArrayXd z(1); z = 1.0;
bool is_propane = (RPname == "PROPANE");
double Tt = (is_propane) ? 85.525 : 273.16,
Tc = (is_propane) ? 369.89 : 647.096;
int NT = 200;
nlohmann::json outputs = nlohmann::json::array();
for (double T : Eigen::ArrayXd::LinSpaced(NT, Tt, Tc)) {
auto o = REFPROP_sat(T);
// Pressure for each phase via REFPROP
double pL = -1; PRESSdll(o.T, o.rhoLmol_L, &(z[0]), pL);
double RL = -1; RMIX2dll(&(z[0]), RL);
double ZLRP = pL/(o.rhoLmol_L*RL*o.T); // Units cancel (factor of 1000 in pL and RL)
double Tr = -1, Dr = -1;
REDXdll(&(z[0]), Tr, Dr);
double pV = -1; PRESSdll(o.T, o.rhoVmol_L, &(z[0]), pV);
double RV = -1; RMIX2dll(&(z[0]), RV);
double ZVRP = pV/(o.rhoVmol_L*RV*o.T); // Units cancel (factor of 1000 in pV and RV)
int itau = 0, idelta = 3; double tau = Tr / o.T, deltaL = o.rhoLmol_L / Dr, Ar03LRP = -1;
double deltaV = o.rhoVmol_L / Dr, Ar03VRP = -1;
PHIXdll(itau, idelta, tau, deltaL, &(z[0]), Ar03LRP);
PHIXdll(itau, idelta, tau, deltaV, &(z[0]), Ar03VRP);
double Ar01LRP = -1, Ar01VRP = -1;
idelta = 1;
PHIXdll(itau, idelta, tau, deltaL, &(z[0]), Ar01LRP);
PHIXdll(itau, idelta, tau, deltaV, &(z[0]), Ar01VRP);
double Ar02LRP = -1, Ar02VRP = -1;
idelta = 2;
PHIXdll(itau, idelta, tau, deltaL, &(z[0]), Ar02LRP);
PHIXdll(itau, idelta, tau, deltaV, &(z[0]), Ar02VRP);
double rhoL = o.rhoLmol_L * 1000.0, rhoV = o.rhoVmol_L*1000.0; for (double Q : { 0, 1 }) {
double rho = (Q == 0) ? rhoL : rhoV;
auto c = with_teqp_and_boost(model, T, rho, z, is_propane);
double Zratiominus1 = c.Zteqp / c.Zexact - 1, Ar01ratiominus1 = c.Ar01teqp / c.Ar01exact - 1;
outputs.push_back({
{"T / K", T},
{"Q", Q},
{"rho / mol/m^3", rho},
{"Zteqp", c.Zteqp},
{"Zexact", c.Zexact},
{"ratio-1", Zratiominus1},
{"ZRP", ((Q == 0) ? ZLRP : ZVRP)},
{"Ar03teqp", c.Ar03teqp},
{"Ar03exact", c.Ar03exact},
{"ratio03-1", c.Ar03teqp / c.Ar03exact - 1},
{"Ar01teqp", c.Ar01teqp},
{"Ar01exact", c.Ar01exact},
{"ratio01-1", Ar01ratiominus1},
{"Ar02teqp", c.Ar02teqp},
{"Ar02exact", c.Ar02exact},
{"ratio02-1", Ar02ratiominus1},
{"Ar01RP", ((Q == 0) ? Ar01LRP : Ar01VRP)},
{"Ar02RP", ((Q == 0) ? Ar02LRP : Ar02VRP)},
{"Ar03RP", ((Q == 0) ? Ar03LRP : Ar03VRP)},
});
}
std::cout << "Completed:" << T << std::endl;
}
std::ofstream file(RPname + "_saturation_Z_accuracy.json");
file << outputs;
return EXIT_SUCCESS;
}
int main()
{
do_one("PROPANE", "n-Propane");
do_one("WATER", "Water");
return EXIT_SUCCESS;
}