2 #if !defined(NO_TABULAR_BACKENDS)
13 static const double Ainv_data[16 * 16] = {
14 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, -2,
15 -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, -2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
16 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, -2, -1, 0,
17 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, -2, 0, 0, 1, 1, 0, 0, -3, 0, 3, 0, 0, 0, 0, 0, -2, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
18 -3, 0, 3, 0, 0, 0, 0, 0, -2, 0, -1, 0, 9, -9, -9, 9, 6, 3, -6, -3, 6, -6, 3, -3, 4, 2, 2, 1, -6, 6, 6, -6, -3, -3, 3, 3, -4,
19 4, -2, 2, -2, -2, -1, -1, 2, 0, -2, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, -2, 0, 0, 0, 0, 0, 1, 0,
20 1, 0, -6, 6, 6, -6, -4, -2, 4, 2, -3, 3, -3, 3, -2, -1, -2, -1, 4, -4, -4, 4, 2, 2, -2, -2, 2, -2, 2, -2, 1, 1, 1, 1};
21 static Eigen::Matrix<double, 16, 16> Ainv(Ainv_data);
34 void load_table(T& table,
const std::string& path_to_tables,
const std::string& filename) {
37 std::string path_to_table = path_to_tables +
"/" + filename;
39 std::cout <<
format(
"Loading table: %s", path_to_table.c_str()) << std::endl;
41 std::vector<char> raw;
45 std::string err =
format(
"Unable to load file %s", path_to_table.c_str());
47 std::cout <<
"err:" << err << std::endl;
51 std::vector<unsigned char> newBuffer(raw.size() * 5);
52 uLong newBufferSize =
static_cast<uLong>(newBuffer.size());
56 code =
uncompress((
unsigned char*)(&(newBuffer[0])), &newBufferSize, (
unsigned char*)(&(raw[0])), rawBufferSize);
59 newBuffer.resize(newBuffer.size() * 5);
60 newBufferSize =
static_cast<uLong>(newBuffer.size());
61 }
else if (code != 0) {
62 std::string err =
format(
"Unable to uncompress file %s with miniz code %d", path_to_table.c_str(), code);
64 std::cout <<
"uncompress err:" << err << std::endl;
70 std::vector<char> charbuffer(newBuffer.begin(), newBuffer.begin() + newBufferSize);
72 msgpack::unpacked msg;
73 msgpack::unpack(msg, &(charbuffer[0]), charbuffer.size());
74 msgpack::object deserialized = msg.get();
77 table.deserialize(deserialized);
80 std::cout <<
format(
"Loaded table: %s in %g sec.", path_to_table.c_str(), (toc - tic) / CLOCKS_PER_SEC) << std::endl;
82 }
catch (std::exception& e) {
83 std::string err =
format(
"Unable to msgpack deserialize %s; err: %s", path_to_table.c_str(), e.what());
85 std::cout <<
"err: " << err << std::endl;
91 void write_table(
const T& table,
const std::string& path_to_tables,
const std::string& name) {
92 msgpack::sbuffer sbuf;
93 msgpack::pack(sbuf, table);
94 std::string tabPath = std::string(path_to_tables +
"/" + name +
".bin");
95 std::string zPath = tabPath +
".z";
96 std::vector<char> buffer(sbuf.size());
97 uLong outSize =
static_cast<uLong>(buffer.size());
98 compress((
unsigned char*)(&(buffer[0])), &outSize, (
unsigned char*)(sbuf.data()),
static_cast<mz_ulong>(sbuf.size()));
99 std::ofstream ofs2(zPath.c_str(), std::ofstream::binary);
100 ofs2.write(&buffer[0], outSize);
104 std::ofstream ofs(tabPath.c_str(), std::ofstream::binary);
105 ofs.write(sbuf.data(), sbuf.size());
114 std::cout <<
format(
"***********************************************\n");
115 std::cout <<
format(
" Saturation Table (%s) \n",
AS->name().c_str());
116 std::cout <<
format(
"***********************************************\n");
125 CoolPropDbl p, pmin = p_triple, pmax = 0.9999 *
AS->p_critical();
126 for (std::size_t i = 0; i <
N - 1; ++i) {
131 p = exp(log(pmin) + (log(pmax) - log(pmin)) / (
N - 1) * i);
138 rhomolarL[i] =
AS->rhomolar();
139 hmolarL[i] =
AS->hmolar();
140 smolarL[i] =
AS->smolar();
141 umolarL[i] =
AS->umolar();
143 logrhomolarL[i] = log(rhomolarL[i]);
144 cpmolarL[i] =
AS->cpmolar();
145 cvmolarL[i] =
AS->cvmolar();
146 speed_soundL[i] =
AS->speed_sound();
147 }
catch (std::exception& e) {
150 std::cout <<
" " << e.what() << std::endl;
156 viscL[i] =
AS->viscosity();
157 condL[i] =
AS->conductivity();
158 logviscL[i] = log(viscL[i]);
159 }
catch (std::exception& e) {
161 std::cout <<
" " << e.what() << std::endl;
169 rhomolarV[i] =
AS->rhomolar();
170 hmolarV[i] =
AS->hmolar();
171 smolarV[i] =
AS->smolar();
172 umolarV[i] =
AS->umolar();
174 logrhomolarV[i] = log(rhomolarV[i]);
175 cpmolarV[i] =
AS->cpmolar();
176 cvmolarV[i] =
AS->cvmolar();
177 speed_soundV[i] =
AS->speed_sound();
178 }
catch (std::exception& e) {
181 std::cout <<
" " << e.what() << std::endl;
187 viscV[i] =
AS->viscosity();
188 condV[i] =
AS->conductivity();
189 logviscV[i] = log(viscV[i]);
190 }
catch (std::exception& e) {
192 std::cout <<
" " << e.what() << std::endl;
201 std::size_t i =
N - 1;
205 rhomolarV[i] =
AS->rhomolar();
206 hmolarV[i] =
AS->hmolar();
207 smolarV[i] =
AS->smolar();
208 umolarV[i] =
AS->umolar();
212 rhomolarL[i] =
AS->rhomolar();
213 hmolarL[i] =
AS->hmolar();
214 smolarL[i] =
AS->smolar();
215 umolarL[i] =
AS->umolar();
217 logpV[i] = log(
AS->p());
218 logrhomolarV[i] = log(rhomolarV[i]);
220 logpL[i] = log(
AS->p());
221 logrhomolarL[i] = log(rhomolarL[i]);
231 std::cout <<
format(
"***********************************************\n");
232 std::cout <<
format(
" Single-Phase Table (%s) \n",
strjoin(AS->fluid_names(),
"&").c_str());
233 std::cout <<
format(
"***********************************************\n");
238 for (std::size_t i = 0; i < Nx; ++i) {
242 x = exp(log(xmin) + (log(xmax) - log(xmin)) / (Nx - 1) * i);
245 x = xmin + (xmax - xmin) / (Nx - 1) * i;
248 for (std::size_t j = 0; j < Ny; ++j) {
252 y = exp(log(ymin) + (log(ymax / ymin)) / (Ny - 1) * j);
255 y = ymin + (ymax - ymin) / (Ny - 1) * j;
260 std::cout <<
"x: " << x <<
" y: " << y << std::endl;
271 AS->update(input_pair, v1, v2);
275 }
catch (std::exception& e) {
278 std::cout <<
" " << e.what() << std::endl;
286 std::cout <<
" 2Phase" << std::endl;
296 rhomolar[i][j] = AS->rhomolar();
297 hmolar[i][j] = AS->hmolar();
298 smolar[i][j] = AS->smolar();
299 umolar[i][j] = AS->umolar();
305 visc[i][j] = AS->viscosity();
306 cond[i][j] = AS->conductivity();
307 }
catch (std::exception&) {
314 dTdx[i][j] = AS->first_partial_deriv(
iT, xkey, ykey);
315 dTdy[i][j] = AS->first_partial_deriv(
iT, ykey, xkey);
316 dpdx[i][j] = AS->first_partial_deriv(
iP, xkey, ykey);
317 dpdy[i][j] = AS->first_partial_deriv(
iP, ykey, xkey);
318 drhomolardx[i][j] = AS->first_partial_deriv(
iDmolar, xkey, ykey);
319 drhomolardy[i][j] = AS->first_partial_deriv(
iDmolar, ykey, xkey);
320 dhmolardx[i][j] = AS->first_partial_deriv(
iHmolar, xkey, ykey);
321 dhmolardy[i][j] = AS->first_partial_deriv(
iHmolar, ykey, xkey);
322 dsmolardx[i][j] = AS->first_partial_deriv(
iSmolar, xkey, ykey);
323 dsmolardy[i][j] = AS->first_partial_deriv(
iSmolar, ykey, xkey);
324 dumolardx[i][j] = AS->first_partial_deriv(
iUmolar, xkey, ykey);
325 dumolardy[i][j] = AS->first_partial_deriv(
iUmolar, ykey, xkey);
330 d2Tdx2[i][j] = AS->second_partial_deriv(
iT, xkey, ykey, xkey, ykey);
331 d2Tdxdy[i][j] = AS->second_partial_deriv(
iT, xkey, ykey, ykey, xkey);
332 d2Tdy2[i][j] = AS->second_partial_deriv(
iT, ykey, xkey, ykey, xkey);
333 d2pdx2[i][j] = AS->second_partial_deriv(
iP, xkey, ykey, xkey, ykey);
334 d2pdxdy[i][j] = AS->second_partial_deriv(
iP, xkey, ykey, ykey, xkey);
335 d2pdy2[i][j] = AS->second_partial_deriv(
iP, ykey, xkey, ykey, xkey);
336 d2rhomolardx2[i][j] = AS->second_partial_deriv(
iDmolar, xkey, ykey, xkey, ykey);
337 d2rhomolardxdy[i][j] = AS->second_partial_deriv(
iDmolar, xkey, ykey, ykey, xkey);
338 d2rhomolardy2[i][j] = AS->second_partial_deriv(
iDmolar, ykey, xkey, ykey, xkey);
339 d2hmolardx2[i][j] = AS->second_partial_deriv(
iHmolar, xkey, ykey, xkey, ykey);
340 d2hmolardxdy[i][j] = AS->second_partial_deriv(
iHmolar, xkey, ykey, ykey, xkey);
341 d2hmolardy2[i][j] = AS->second_partial_deriv(
iHmolar, ykey, xkey, ykey, xkey);
342 d2smolardx2[i][j] = AS->second_partial_deriv(
iSmolar, xkey, ykey, xkey, ykey);
343 d2smolardxdy[i][j] = AS->second_partial_deriv(
iSmolar, xkey, ykey, ykey, xkey);
344 d2smolardy2[i][j] = AS->second_partial_deriv(
iSmolar, ykey, xkey, ykey, xkey);
345 d2umolardx2[i][j] = AS->second_partial_deriv(
iUmolar, xkey, ykey, xkey, ykey);
346 d2umolardxdy[i][j] = AS->second_partial_deriv(
iUmolar, xkey, ykey, ykey, xkey);
347 d2umolardy2[i][j] = AS->second_partial_deriv(
iUmolar, ykey, xkey, ykey, xkey);
352 std::vector<std::string> fluids = AS->fluid_names();
353 std::vector<CoolPropDbl> fractions = AS->get_mole_fractions();
354 std::vector<std::string> components;
355 for (std::size_t i = 0; i < fluids.size(); ++i) {
356 components.push_back(
format(
"%s[%0.10Lf]", fluids[i].c_str(), fractions[i]));
358 std::string table_directory =
get_home_dir() +
"/.CoolProp/Tables/";
359 std::string alt_table_directory =
get_config_string(ALTERNATIVE_TABLES_DIRECTORY);
360 if (!alt_table_directory.empty()) {
361 table_directory = alt_table_directory;
363 return table_directory + AS->backend_name() +
"(" +
strjoin(components,
"&") +
")";
367 std::string path_to_tables = this->path_to_tables();
375 write_table(single_phase_logph, path_to_tables,
"single_phase_logph");
376 write_table(single_phase_logpT, path_to_tables,
"single_phase_logpT");
377 write_table(pure_saturation, path_to_tables,
"pure_saturation");
378 write_table(phase_envelope, path_to_tables,
"phase_envelope");
383 if (loaded ==
false) {
387 std::cout <<
"Tables loaded" << std::endl;
397 return phase_envelope_sat(phase_envelope, key,
iP, _p) * factor;
399 return pure_saturation.
evaluate(key, _p, 1, cached_saturation_iL, cached_saturation_iV) * factor;
408 return phase_envelope_sat(phase_envelope, key,
iP, _p) * factor;
410 return pure_saturation.
evaluate(key, _p, 0, cached_saturation_iL, cached_saturation_iV) * factor;
416 if (using_single_phase_table) {
420 return phase_envelope_sat(phase_envelope,
iP,
iT, _T);
429 if (using_single_phase_table) {
430 switch (selected_table) {
431 case SELECTED_PH_TABLE:
432 return evaluate_single_phase_phmolar(
iT, cached_single_phase_i, cached_single_phase_j);
433 case SELECTED_PT_TABLE:
435 case SELECTED_NO_TABLE:
441 return phase_envelope_sat(phase_envelope,
iT,
iP, _p);
446 return pure_saturation.
evaluate(
iT, _p, _Q, cached_saturation_iL, cached_saturation_iV);
454 if (using_single_phase_table) {
455 switch (selected_table) {
456 case SELECTED_PH_TABLE:
457 return evaluate_single_phase_phmolar(
iDmolar, cached_single_phase_i, cached_single_phase_j);
458 case SELECTED_PT_TABLE:
459 return evaluate_single_phase_pT(
iDmolar, cached_single_phase_i, cached_single_phase_j);
460 case SELECTED_NO_TABLE:
466 return phase_envelope_sat(phase_envelope,
iDmolar,
iP, _p);
468 return pure_saturation.
evaluate(
iDmolar, _p, _Q, cached_saturation_iL, cached_saturation_iV);
475 if (using_single_phase_table) {
476 switch (selected_table) {
477 case SELECTED_PH_TABLE:
479 case SELECTED_PT_TABLE:
480 return evaluate_single_phase_pT(
iHmolar, cached_single_phase_i, cached_single_phase_j);
481 case SELECTED_NO_TABLE:
487 return phase_envelope_sat(phase_envelope,
iHmolar,
iP, _p);
489 return pure_saturation.
evaluate(
iHmolar, _p, _Q, cached_saturation_iL, cached_saturation_iV);
496 if (using_single_phase_table) {
497 switch (selected_table) {
498 case SELECTED_PH_TABLE:
499 return evaluate_single_phase_phmolar(
iSmolar, cached_single_phase_i, cached_single_phase_j);
500 case SELECTED_PT_TABLE:
501 return evaluate_single_phase_pT(
iSmolar, cached_single_phase_i, cached_single_phase_j);
502 case SELECTED_NO_TABLE:
508 return phase_envelope_sat(phase_envelope,
iSmolar,
iP, _p);
510 return pure_saturation.
evaluate(
iSmolar, _p, _Q, cached_saturation_iL, cached_saturation_iV);
517 if (using_single_phase_table) {
518 switch (selected_table) {
519 case SELECTED_PH_TABLE:
520 return evaluate_single_phase_phmolar(
iUmolar, cached_single_phase_i, cached_single_phase_j);
521 case SELECTED_PT_TABLE:
522 return evaluate_single_phase_pT(
iUmolar, cached_single_phase_i, cached_single_phase_j);
523 case SELECTED_NO_TABLE:
532 return hmolar - _p / rhomolar;
534 return pure_saturation.
evaluate(
iUmolar, _p, _Q, cached_saturation_iL, cached_saturation_iV);
541 if (using_single_phase_table) {
545 return phase_envelope_sat(phase_envelope,
iCpmolar,
iP, _p);
547 return pure_saturation.
evaluate(
iCpmolar, _p, _Q, cached_saturation_iL, cached_saturation_iV);
554 if (using_single_phase_table) {
558 return phase_envelope_sat(phase_envelope,
iCvmolar,
iP, _p);
560 return pure_saturation.
evaluate(
iCvmolar, _p, _Q, cached_saturation_iL, cached_saturation_iV);
568 if (using_single_phase_table) {
569 switch (selected_table) {
570 case SELECTED_PH_TABLE:
571 return evaluate_single_phase_phmolar_transport(
iviscosity, cached_single_phase_i, cached_single_phase_j);
572 case SELECTED_PT_TABLE:
573 return evaluate_single_phase_pT_transport(
iviscosity, cached_single_phase_i, cached_single_phase_j);
574 case SELECTED_NO_TABLE:
580 return phase_envelope_sat(phase_envelope,
iviscosity,
iP, _p);
582 return pure_saturation.
evaluate(
iviscosity, _p, _Q, cached_saturation_iL, cached_saturation_iV);
589 if (using_single_phase_table) {
590 switch (selected_table) {
591 case SELECTED_PH_TABLE:
592 return evaluate_single_phase_phmolar_transport(
iconductivity, cached_single_phase_i, cached_single_phase_j);
593 case SELECTED_PT_TABLE:
594 return evaluate_single_phase_pT_transport(
iconductivity, cached_single_phase_i, cached_single_phase_j);
595 case SELECTED_NO_TABLE:
610 if (using_single_phase_table) {
611 return sqrt(1 / molar_mass() * cpmolar() / cvmolar() * first_partial_deriv(
iP,
iDmolar,
iT));
616 return pure_saturation.
evaluate(
ispeed_sound, _p, _Q, cached_saturation_iL, cached_saturation_iV);
621 if (using_single_phase_table) {
622 CoolPropDbl dOf_dx, dOf_dy, dWrt_dx, dWrt_dy, dConstant_dx, dConstant_dy;
625 double Of_conversion_factor = 1.0, Wrt_conversion_factor = 1.0, Constant_conversion_factor = 1.0, MM = AS->molar_mass();
630 switch (selected_table) {
631 case SELECTED_PH_TABLE: {
632 dOf_dx = evaluate_single_phase_phmolar_derivative(Of, cached_single_phase_i, cached_single_phase_j, 1, 0);
633 dOf_dy = evaluate_single_phase_phmolar_derivative(Of, cached_single_phase_i, cached_single_phase_j, 0, 1);
634 dWrt_dx = evaluate_single_phase_phmolar_derivative(Wrt, cached_single_phase_i, cached_single_phase_j, 1, 0);
635 dWrt_dy = evaluate_single_phase_phmolar_derivative(Wrt, cached_single_phase_i, cached_single_phase_j, 0, 1);
636 dConstant_dx = evaluate_single_phase_phmolar_derivative(Constant, cached_single_phase_i, cached_single_phase_j, 1, 0);
637 dConstant_dy = evaluate_single_phase_phmolar_derivative(Constant, cached_single_phase_i, cached_single_phase_j, 0, 1);
640 case SELECTED_PT_TABLE: {
641 dOf_dx = evaluate_single_phase_pT_derivative(Of, cached_single_phase_i, cached_single_phase_j, 1, 0);
642 dOf_dy = evaluate_single_phase_pT_derivative(Of, cached_single_phase_i, cached_single_phase_j, 0, 1);
643 dWrt_dx = evaluate_single_phase_pT_derivative(Wrt, cached_single_phase_i, cached_single_phase_j, 1, 0);
644 dWrt_dy = evaluate_single_phase_pT_derivative(Wrt, cached_single_phase_i, cached_single_phase_j, 0, 1);
645 dConstant_dx = evaluate_single_phase_pT_derivative(Constant, cached_single_phase_i, cached_single_phase_j, 1, 0);
646 dConstant_dy = evaluate_single_phase_pT_derivative(Constant, cached_single_phase_i, cached_single_phase_j, 0, 1);
649 case SELECTED_NO_TABLE:
652 double val = (dOf_dx * dConstant_dy - dOf_dy * dConstant_dx) / (dWrt_dx * dConstant_dy - dWrt_dy * dConstant_dx);
653 return val * Of_conversion_factor / Wrt_conversion_factor;
655 throw ValueError(
format(
"Inputs [rho: %g mol/m3, T: %g K, p: %g Pa] are two-phase; cannot use single-phase derivatives", _rhomolar, _T, _p));
661 if (AS->get_mole_fractions().size() > 1) {
662 throw ValueError(
"calc_first_saturation_deriv not available for mixtures");
664 if (std::abs(_Q) < 1e-6) {
666 }
else if (std::abs(_Q - 1) < 1e-6) {
669 throw ValueError(
format(
"Quality [%Lg] must be either 0 or 1 to within 1 ppm", _Q));
679 return -
POW2(rhomolar()) * (1 / rhoV - 1 / rhoL) / (hV - hL);
694 CoolPropDbl dxdp_h = (Q() * dhV_dp + (1 - Q()) * dhL_dp) / (hL - hV);
695 CoolPropDbl dvdp_h = dvL_dp + dxdp_h * (1 / rhoV - 1 / rhoL) + Q() * (dvV_dp - dvL_dp);
696 return -
POW2(rhomolar()) * dvdp_h;
700 throw ValueError(
"These inputs are not supported to calc_first_two_phase_deriv");
708 bool drho_dh__p =
false;
709 bool drho_dp__h =
false;
710 bool rho_spline =
false;
714 if (_drho_spline_dh__constp)
return _drho_spline_dh__constp;
719 if (_drho_spline_dp__consth)
return _drho_spline_dp__consth;
721 return first_two_phase_deriv_splined(
iDmolar,
iP,
iHmolar, x_end) * molar_mass();
726 if (_rho_spline)
return _rho_spline;
730 throw ValueError(
"These inputs are not supported to calc_first_two_phase_deriv");
734 throw ValueError(
format(
"Q [%g] is greater than x_end [%Lg]", _Q, x_end).c_str());
764 CoolPropDbl TE = pure_saturation.
evaluate(
iT, _p, x_end, cached_saturation_iL, cached_saturation_iV);
779 CoolPropDbl Abracket = (2 * dL - 2 * dE + Delta_end * (drho_dh_liq__constp + drho_dh_end__constp));
781 CoolPropDbl b = 3 /
POW2(Delta_end) * (-dL + dE) - 1 / Delta_end * (drho_dh_end__constp + 2 * drho_dh_liq__constp);
786 _rho_spline = a *
POW3(Delta) + b *
POW2(Delta) + c * Delta + d;
787 _drho_spline_dh__constp = 3 * a *
POW2(Delta) + 2 * b * Delta + c;
788 if (rho_spline)
return _rho_spline;
789 if (drho_dh__p)
return _drho_spline_dh__constp;
801 CoolPropDbl drho_dp_end =
POW2(dE) * (x_end /
POW2(dV) * drhoV_dp_sat + (1 - x_end) /
POW2(dL) * drhoL_dp_sat);
815 CoolPropDbl d_Delta_end_dp__consth = x_end * (dhV_dp_sat - dhL_dp_sat);
819 CoolPropDbl d_Abracket_dp_consth = (2 * drhoL_dp_sat - 2 * drho_dp_end + Delta_end * (d2rhodhdp_liq + d2rhodhdp_end)
820 + d_Delta_end_dp__consth * (drho_dh_liq__constp + drho_dh_end__constp));
821 CoolPropDbl da_dp = 1 /
POW3(Delta_end) * d_Abracket_dp_consth + Abracket * (-3 /
POW4(Delta_end) * d_Delta_end_dp__consth);
822 CoolPropDbl db_dp = -6 /
POW3(Delta_end) * d_Delta_end_dp__consth * (dE - dL) + 3 /
POW2(Delta_end) * (drho_dp_end - drhoL_dp_sat)
823 + (1 /
POW2(Delta_end) * d_Delta_end_dp__consth) * (drho_dh_end__constp + 2 * drho_dh_liq__constp)
824 - (1 / Delta_end) * (d2rhodhdp_end + 2 * d2rhodhdp_liq);
828 _drho_spline_dp__consth =
829 (3 * a *
POW2(Delta) + 2 * b * Delta + c) * d_Delta_dp__consth +
POW3(Delta) * da_dp +
POW2(Delta) * db_dp + Delta * dc_dp + dd_dp;
830 if (drho_dp__h)
return _drho_spline_dp__consth;
832 throw ValueError(
"Something went wrong in TabularBackend::calc_first_two_phase_deriv_splined");
847 mass_to_molar_inputs(input_pair, ld_value1, ld_value2);
855 cached_single_phase_i = std::numeric_limits<std::size_t>::max();
856 cached_single_phase_j = std::numeric_limits<std::size_t>::max();
857 cached_saturation_iL = std::numeric_limits<std::size_t>::max();
858 cached_saturation_iV = std::numeric_limits<std::size_t>::max();
868 switch (input_pair) {
874 using_single_phase_table =
false;
876 std::cout <<
"inputs are not in range";
880 using_single_phase_table =
true;
881 std::size_t iL = std::numeric_limits<std::size_t>::max(), iV = std::numeric_limits<std::size_t>::max(), iclosest = 0;
884 bool is_two_phase =
false;
897 using_single_phase_table =
false;
898 _Q = (
static_cast<double>(_hmolar) - hL) / (hV - hL);
901 format(
"vapor quality is not in (0,1) for hmolar: %g p: %g, hL: %g hV: %g ",
static_cast<double>(_hmolar), _p, hL, hV));
903 cached_saturation_iL = iL;
904 cached_saturation_iV = iV;
908 selected_table = SELECTED_PH_TABLE;
910 find_native_nearest_good_indices(single_phase_logph, dataset->coeffs_ph, _hmolar, _p, cached_single_phase_i,
911 cached_single_phase_j);
913 recalculate_singlephase_phase();
923 using_single_phase_table =
false;
925 std::cout <<
"inputs are not in range";
927 throw ValueError(
format(
"inputs are not in range, p=%g Pa, T=%g K", _p, _T));
929 using_single_phase_table =
true;
930 std::size_t iL = 0, iV = 0, iclosest = 0;
933 bool is_two_phase =
false;
941 is_two_phase = pure_saturation.
is_inside(
iP, _p,
iT, _T, iL, iV, TL, TV);
945 using_single_phase_table =
false;
948 selected_table = SELECTED_PT_TABLE;
950 find_native_nearest_good_indices(single_phase_logpT, dataset->coeffs_pT, _T, _p, cached_single_phase_i, cached_single_phase_j);
953 double rhoc = rhomolar_critical();
954 if (imposed_phase_index ==
iphase_liquid && cached_single_phase_i > 0) {
956 if (_p < this->AS->p_critical()) {
957 while (cached_single_phase_i > 0
958 && single_phase_logpT.rhomolar[cached_single_phase_i + 1][cached_single_phase_j] < rhoc) {
960 cached_single_phase_i--;
962 double rho = evaluate_single_phase_pT(
iDmolar, cached_single_phase_i, cached_single_phase_j);
971 && cached_single_phase_i > 0) {
974 if (_p < this->AS->p_critical()) {
975 while (cached_single_phase_i > 0
976 && single_phase_logpT.rhomolar[cached_single_phase_i][cached_single_phase_j + 1] > rhoc) {
978 cached_single_phase_i++;
980 double rho = evaluate_single_phase_pT(
iDmolar, cached_single_phase_i, cached_single_phase_j);
994 if (_p < this->AS->p_critical()) {
995 double Ts = pure_saturation.
evaluate(
iT, _p, _Q, iL, iV);
996 double TL = single_phase_logpT.T[cached_single_phase_i][cached_single_phase_j];
997 double TR = single_phase_logpT.T[cached_single_phase_i + 1][cached_single_phase_j];
998 if (TL < Ts && Ts < TR) {
1000 if (cached_single_phase_i == 0) {
1001 throw ValueError(
format(
"P, T are near saturation, but cannot move the cell to the left"));
1004 cached_single_phase_i--;
1006 if (cached_single_phase_i > single_phase_logpT.
Nx - 2) {
1007 throw ValueError(
format(
"P,T are near saturation, but cannot move the cell to the right"));
1010 cached_single_phase_i++;
1016 recalculate_singlephase_phase();
1026 switch (input_pair) {
1049 using_single_phase_table =
true;
1050 std::size_t iL = std::numeric_limits<std::size_t>::max(), iV = std::numeric_limits<std::size_t>::max();
1052 std::size_t iclosest = 0;
1054 bool is_two_phase =
false;
1062 std::vector<std::pair<std::size_t, std::size_t>> intersect =
1064 if (intersect.size() < 2) {
1067 iV = intersect[0].first;
1068 iL = intersect[1].first;
1073 is_two_phase = pure_saturation.
is_inside(
iP, _p, otherkey, otherval, iL, iV, zL, zV);
1078 using_single_phase_table =
false;
1080 _Q = (1 / otherval - 1 / zL) / (1 / zV - 1 / zL);
1082 _Q = (otherval - zL) / (zV - zL);
1087 otherval,
static_cast<double>(_p)));
1088 }
else if (!is_mixture) {
1089 cached_saturation_iL = iL;
1090 cached_saturation_iV = iV;
1094 selected_table = SELECTED_PH_TABLE;
1096 find_nearest_neighbor(single_phase_logph, dataset->coeffs_ph,
iP, _p, otherkey, otherval, cached_single_phase_i,
1097 cached_single_phase_j);
1099 invert_single_phase_x(single_phase_logph, dataset->coeffs_ph, otherkey, otherval, _p, cached_single_phase_i, cached_single_phase_j);
1101 recalculate_singlephase_phase();
1109 switch (input_pair) {
1126 using_single_phase_table =
true;
1127 std::size_t iL = std::numeric_limits<std::size_t>::max(), iV = std::numeric_limits<std::size_t>::max();
1129 std::size_t iclosest = 0;
1131 bool is_two_phase =
false;
1139 std::vector<std::pair<std::size_t, std::size_t>> intersect =
1141 if (intersect.size() < 2) {
1144 iV = intersect[0].first;
1145 iL = intersect[1].first;
1150 is_two_phase = pure_saturation.
is_inside(
iT, _T, otherkey, otherval, iL, iV, zL, zV);
1154 using_single_phase_table =
false;
1156 _Q = (1 / otherval - 1 / zL) / (1 / zV - 1 / zL);
1158 _Q = (otherval - zL) / (zV - zL);
1162 otherval,
static_cast<double>(_T)));
1163 }
else if (!is_mixture) {
1164 cached_saturation_iL = iL;
1165 cached_saturation_iV = iV;
1166 _p = pure_saturation.
evaluate(
iP, _T, _Q, iL, iV);
1170 if (intersect.size() < 2) {
1173 iV = intersect[0].first;
1174 iL = intersect[1].first;
1177 _p = _Q * pV + (1 - _Q) * pL;
1180 selected_table = SELECTED_PT_TABLE;
1182 find_nearest_neighbor(single_phase_logpT, dataset->coeffs_pT,
iT, _T, otherkey, otherval, cached_single_phase_i,
1183 cached_single_phase_j);
1185 invert_single_phase_y(single_phase_logpT, dataset->coeffs_pT, otherkey, otherval, _T, cached_single_phase_i, cached_single_phase_j);
1187 recalculate_singlephase_phase();
1192 std::size_t iL = 0, iV = 0;
1195 using_single_phase_table =
false;
1197 throw ValueError(
format(
"vapor quality [%g] is not in (0,1)",
static_cast<double>(val2)));
1202 if (intersect.size() < 2) {
1205 iV = intersect[0].first;
1206 iL = intersect[1].first;
1210 bool it_is_inside = pure_saturation.
is_inside(
iP, _p,
iQ, _Q, iL, iV, TL, TV);
1211 if (!it_is_inside) {
1212 throw ValueError(
"Not possible to determine whether pressure is inside or not");
1215 _T = _Q * TV + (1 - _Q) * TL;
1216 cached_saturation_iL = iL;
1217 cached_saturation_iV = iV;
1223 std::size_t iL = 0, iV = 0;
1227 using_single_phase_table =
false;
1229 throw ValueError(
format(
"vapor quality [%g] is not in (0,1)",
static_cast<double>(val1)));
1234 if (intersect.size() < 2) {
1237 iV = intersect[0].first;
1238 iL = intersect[1].first;
1244 _p = _Q * pV + (1 - _Q) * pL;
1245 cached_saturation_iL = iL;
1246 cached_saturation_iV = iV;
1252 throw ValueError(
"Sorry, but this set of inputs is not supported for Tabular backend");
1258 write_table(single_phase_logph, path_to_tables,
"single_phase_logph");
1259 write_table(single_phase_logpT, path_to_tables,
"single_phase_logpT");
1260 write_table(pure_saturation, path_to_tables,
"pure_saturation");
1261 write_table(phase_envelope, path_to_tables,
"phase_envelope");
1265 single_phase_logph.AS = AS;
1266 single_phase_logpT.AS = AS;
1267 pure_saturation.AS = AS;
1268 single_phase_logph.set_limits();
1269 single_phase_logpT.set_limits();
1270 load_table(single_phase_logph, path_to_tables,
"single_phase_logph.bin.z");
1271 load_table(single_phase_logpT, path_to_tables,
"single_phase_logpT.bin.z");
1272 load_table(pure_saturation, path_to_tables,
"pure_saturation.bin.z");
1273 load_table(phase_envelope, path_to_tables,
"phase_envelope.bin.z");
1274 tables_loaded =
true;
1276 std::cout <<
"Tables loaded" << std::endl;
1282 if (AS->get_mole_fractions().size() == 1) {
1283 pure_saturation.build(AS);
1286 AS->build_phase_envelope(
"");
1290 phase_envelope.copy_from_nonpackable(PED);
1292 pure_saturation.
resize(pure_saturation.N);
1294 single_phase_logph.build(AS);
1295 single_phase_logpT.build(AS);
1296 tables_loaded =
true;
1301 const std::string path = path_to_tables(AS);
1303 std::map<std::string, TabularDataSet>::iterator it = data.find(path);
1305 if (it != data.end()) {
1306 loaded = it->second.tables_loaded;
1307 return &(it->second);
1312 data.insert(std::pair<std::string, TabularDataSet>(path, set));
1319 }
catch (std::exception&) {
1327 if (!coeffs.empty()) {
1331 const int param_count = 6;
1333 std::vector<std::vector<double>>*f = NULL, *fx = NULL, *fy = NULL, *fxy = NULL;
1335 clock_t t1 = clock();
1338 coeffs.resize(table.
Nx - 1, std::vector<CellCoeffs>(table.
Ny - 1));
1340 int valid_cell_count = 0;
1341 for (std::size_t k = 0; k < param_count; ++k) {
1343 if (param == table.
xkey || param == table.
ykey) {
1352 fxy = &(table.d2Tdxdy);
1358 fxy = &(table.d2pdxdy);
1361 f = &(table.rhomolar);
1362 fx = &(table.drhomolardx);
1363 fy = &(table.drhomolardy);
1364 fxy = &(table.d2rhomolardxdy);
1367 f = &(table.smolar);
1368 fx = &(table.dsmolardx);
1369 fy = &(table.dsmolardy);
1370 fxy = &(table.d2smolardxdy);
1373 f = &(table.hmolar);
1374 fx = &(table.dhmolardx);
1375 fy = &(table.dhmolardy);
1376 fxy = &(table.d2hmolardxdy);
1379 f = &(table.umolar);
1380 fx = &(table.dumolardx);
1381 fy = &(table.dumolardy);
1382 fxy = &(table.d2umolardxdy);
1385 throw ValueError(
"Invalid variable type to build_coeffs");
1387 for (std::size_t i = 0; i < table.
Nx - 1; ++i)
1389 for (std::size_t j = 0; j < table.
Ny - 1; ++j)
1394 Eigen::Matrix<double, 16, 1> F;
1397 F(1) = (*f)[i + 1][j];
1398 F(2) = (*f)[i][j + 1];
1399 F(3) = (*f)[i + 1][j + 1];
1402 coeffs[i][j].dx_dxhat = table.
xvec[i + 1] - table.
xvec[i];
1403 double dx_dxhat = coeffs[i][j].dx_dxhat;
1404 F(4) = (*fx)[i][j] * dx_dxhat;
1405 F(5) = (*fx)[i + 1][j] * dx_dxhat;
1406 F(6) = (*fx)[i][j + 1] * dx_dxhat;
1407 F(7) = (*fx)[i + 1][j + 1] * dx_dxhat;
1410 coeffs[i][j].dy_dyhat = table.
yvec[j + 1] - table.
yvec[j];
1411 double dy_dyhat = coeffs[i][j].dy_dyhat;
1412 F(8) = (*fy)[i][j] * dy_dyhat;
1413 F(9) = (*fy)[i + 1][j] * dy_dyhat;
1414 F(10) = (*fy)[i][j + 1] * dy_dyhat;
1415 F(11) = (*fy)[i + 1][j + 1] * dy_dyhat;
1417 F(12) = (*fxy)[i][j] * dy_dyhat * dx_dxhat;
1418 F(13) = (*fxy)[i + 1][j] * dy_dyhat * dx_dxhat;
1419 F(14) = (*fxy)[i][j + 1] * dy_dyhat * dx_dxhat;
1420 F(15) = (*fxy)[i + 1][j + 1] * dy_dyhat * dx_dxhat;
1422 Eigen::MatrixXd alpha = Ainv.transpose() * F;
1424 coeffs[i][j].set(param, valpha);
1425 coeffs[i][j].set_valid();
1428 coeffs[i][j].set_invalid();
1432 double elapsed = (clock() - t1) / ((
double)CLOCKS_PER_SEC);
1434 std::cout <<
format(
"Calculated bicubic coefficients for %d good cells in %g sec.\n", valid_cell_count, elapsed);
1436 std::size_t remap_count = 0;
1438 for (std::size_t i = 0; i < table.
Nx - 1; ++i)
1440 for (std::size_t j = 0; j < table.
Ny - 1; ++j)
1443 if (!coeffs[i][j].valid()) {
1445 int xoffsets[] = {-1, 1, 0, 0, -1, 1, 1, -1};
1446 int yoffsets[] = {0, 0, 1, -1, -1, -1, 1, 1};
1448 std::size_t N =
sizeof(xoffsets) /
sizeof(xoffsets[0]);
1449 for (std::size_t k = 0; k < N; ++k) {
1450 std::size_t iplus = i + xoffsets[k];
1451 std::size_t jplus = j + yoffsets[k];
1452 if (0 < iplus && iplus < table.
Nx - 1 && 0 < jplus && jplus < table.
Ny - 1 && coeffs[iplus][jplus].valid()) {
1453 coeffs[i][j].set_alternate(iplus, jplus);
1456 std::cout <<
format(
"Mapping %d,%d to %d,%d\n", i, j, iplus, jplus);
1465 std::cout <<
format(
"Remapped %d cells\n", remap_count);
1470 # if defined(ENABLE_CATCH)
1471 # include <catch2/catch_all.hpp>
1474 static shared_ptr<CoolProp::AbstractState> ASHEOS, ASTTSE, ASBICUBIC;
1477 class TabularFixture
1482 if (ASHEOS.get() == NULL) {
1485 if (ASTTSE.get() == NULL) {
1488 if (ASBICUBIC.get() == NULL) {
1493 TEST_CASE_METHOD(TabularFixture,
"Tests for tabular backends with water",
"[Tabular]") {
1494 SECTION(
"first_saturation_deriv invalid quality") {
1502 SECTION(
"first_saturation_deriv dp/dT") {
1511 CAPTURE(actual_TTSE);
1512 CAPTURE(actual_BICUBIC);
1513 CHECK(std::abs((expected - actual_TTSE) / expected) < 1e-6);
1514 CHECK(std::abs((expected - actual_BICUBIC) / expected) < 1e-6);
1516 SECTION(
"first_saturation_deriv dDmolar/dP") {
1525 CAPTURE(actual_TTSE);
1526 CAPTURE(actual_BICUBIC);
1527 CHECK(std::abs((expected - actual_TTSE) / expected) < 1e-6);
1528 CHECK(std::abs((expected - actual_BICUBIC) / expected) < 1e-6);
1530 SECTION(
"first_saturation_deriv dDmass/dP") {
1539 CAPTURE(actual_TTSE);
1540 CAPTURE(actual_BICUBIC);
1541 CHECK(std::abs((expected - actual_TTSE) / expected) < 1e-6);
1542 CHECK(std::abs((expected - actual_BICUBIC) / expected) < 1e-6);
1544 SECTION(
"first_saturation_deriv dHmass/dP") {
1553 CAPTURE(actual_TTSE);
1554 CAPTURE(actual_BICUBIC);
1555 CHECK(std::abs((expected - actual_TTSE) / expected) < 1e-6);
1556 CHECK(std::abs((expected - actual_BICUBIC) / expected) < 1e-6);
1558 SECTION(
"first_saturation_deriv dHmass/dP w/ QT as inputs") {
1567 CAPTURE(actual_TTSE);
1568 CAPTURE(actual_BICUBIC);
1569 CHECK(std::abs((expected - actual_TTSE) / expected) < 1e-6);
1570 CHECK(std::abs((expected - actual_BICUBIC) / expected) < 1e-6);
1572 SECTION(
"first_two_phase_deriv dDmolar/dP|Hmolar") {
1581 CAPTURE(actual_TTSE);
1582 CAPTURE(actual_BICUBIC);
1583 CHECK(std::abs((expected - actual_TTSE) / expected) < 1e-6);
1584 CHECK(std::abs((expected - actual_BICUBIC) / expected) < 1e-6);
1586 SECTION(
"first_two_phase_deriv dDmass/dP|Hmass") {
1595 CAPTURE(actual_TTSE);
1596 CAPTURE(actual_BICUBIC);
1597 CHECK(std::abs((expected - actual_TTSE) / expected) < 1e-6);
1598 CHECK(std::abs((expected - actual_BICUBIC) / expected) < 1e-6);
1600 SECTION(
"first_two_phase_deriv dDmass/dHmass|P") {
1609 CAPTURE(actual_TTSE);
1610 CAPTURE(actual_BICUBIC);
1611 CHECK(std::abs((expected - actual_TTSE) / expected) < 1e-6);
1612 CHECK(std::abs((expected - actual_BICUBIC) / expected) < 1e-6);
1614 SECTION(
"first_partial_deriv dHmass/dT|P") {
1617 double expected = ASHEOS->cpmass();
1624 CAPTURE(dhdT_BICUBIC);
1625 CHECK(std::abs((expected - dhdT_TTSE) / expected) < 1e-4);
1626 CHECK(std::abs((expected - dhdT_BICUBIC) / expected) < 1e-4);
1628 SECTION(
"first_partial_deriv dHmolar/dT|P") {
1631 double expected = ASHEOS->cpmolar();
1638 CAPTURE(dhdT_BICUBIC);
1639 CHECK(std::abs((expected - dhdT_TTSE) / expected) < 1e-4);
1640 CHECK(std::abs((expected - dhdT_BICUBIC) / expected) < 1e-4);
1642 SECTION(
"check isentropic process") {
1648 double s0 = ASHEOS->smolar();
1650 double expected = ASHEOS->T();
1652 double actual_TTSE = ASTTSE->T();
1654 double actual_BICUBIC = ASBICUBIC->T();
1656 CAPTURE(actual_TTSE);
1657 CAPTURE(actual_BICUBIC);
1658 CHECK(std::abs((expected - actual_TTSE) / expected) < 1e-2);
1659 CHECK(std::abs((expected - actual_BICUBIC) / expected) < 1e-2);
1661 SECTION(
"check D=1 mol/m3, T=500 K inputs") {
1666 double expected = ASHEOS->p();
1668 double actual_TTSE = ASTTSE->p();
1670 double actual_BICUBIC = ASBICUBIC->p();
1672 CAPTURE(actual_TTSE);
1673 CAPTURE(actual_BICUBIC);
1674 CHECK(std::abs((expected - actual_TTSE) / expected) < 1e-3);
1675 CHECK(std::abs((expected - actual_BICUBIC) / expected) < 1e-3);