1 #ifndef PHASEENVELOPE_H
2 #define PHASEENVELOPE_H
16 throw ValueError(
"Mole fractions have not been set yet.");
25 std::vector<CoolPropDbl> Tbp, Qbp;
26 std::vector<std::size_t> Nbp;
40 if (max_sat_T.
rhomolar < crit.rhomolar && max_sat_T.
rhomolar < max_sat_p.rhomolar) {
42 if (max_sat_p.rhomolar < crit.rhomolar) {
46 Tbp.push_back(max_sat_p.T);
47 Tbp.push_back(crit.T);
53 Tbp.push_back(crit.T);
54 Tbp.push_back(max_sat_p.T);
66 for (std::size_t i = 0; i < Tbp.size() - 1; ++i) {
68 std::size_t N = Nbp[i];
79 std::vector<CoolPropDbl>(1, 1.0), std::vector<CoolPropDbl>(1, 1.0));
84 std::vector<CoolPropDbl>(1, 1.0), std::vector<CoolPropDbl>(1, 1.0));
94 std::vector<CriticalState> critpts;
104 std::size_t failure_count = 0;
150 NR.
call(HEOS, IO.
y, IO.
x, IO);
155 bool dont_extrapolate =
false;
160 std::size_t iter = 0,
167 if (failure_count > 5) {
173 if (iter - iter0 > 0) {
176 if (dont_extrapolate) {
179 }
else if (iter - iter0 == 2) {
182 for (std::size_t i = 0; i < IO.
x.size() - 1; ++i)
186 }
else if (iter - iter0 == 3) {
189 for (std::size_t i = 0; i < IO.
x.size() - 1; ++i)
193 }
else if (iter - iter0 > 3) {
202 if (std::abs((T_linear - IO.
T) / IO.
T) > 0.1) {
205 factor = 1 + (factor - 1) / 2;
209 for (std::size_t i = 0; i < IO.
x.size() - 1; ++i)
215 if (IO.
x[i] < 0 || IO.
x[i] > 1) {
218 factor = 1 + (factor - 1) / 2;
226 IO.
x[IO.
x.size() - 1] = 1 - std::accumulate(IO.
x.begin(), IO.
x.end() - 1, 0.0);
233 NR.
call(HEOS, IO.
y, IO.
x, IO);
246 if (!env.T.empty() && std::abs(env.T[env.T.size() - 1] - IO.
T) > 100) {
247 throw ValueError(
"Change in temperature too large");
249 }
catch (std::exception& e) {
251 std::cout << e.what() << std::endl;
257 throw ValueError(
format(
"Unable to calculate at least 4 points in phase envelope; quitting"));
260 factor = 1 + (factor - 1) / 2;
268 <<
" Ns " << IO.
Nsteps <<
" factor " << factor << std::endl;
311 dont_extrapolate =
false;
316 factor = 1 + (factor - 1) / 10;
317 }
else if (IO.
Nsteps > 5) {
318 factor = 1 + (factor - 1) / 3;
319 }
else if (IO.
Nsteps <= 4) {
320 factor = 1 + (factor - 1) * 2;
323 factor = std::max(factor,
static_cast<CoolPropDbl>(1.01));
327 factor = std::min(factor,
static_cast<CoolPropDbl>(1.1));
332 CoolPropDbl max_fraction = *std::max_element(IO.
x.begin(), IO.
x.end());
333 if (iter > 4 && (IO.
p < env.p[0] || std::abs(1.0 - max_fraction) < 1e-9)) {
336 std::cout <<
format(
"envelope built.\n");
337 std::cout <<
format(
"closest fraction to 1.0: distance %g\n", 1 - max_fraction);
361 double acceptable_pdiff = 0.5;
362 double acceptable_rhodiff = 0.25;
364 if (level ==
"veryfine") {
365 acceptable_pdiff = 0.1;
366 acceptable_rhodiff = 0.1;
368 if (level ==
"none") {
375 if ((std::abs(env.rhomolar_vap[i] / env.rhomolar_vap[i + 1] - 1) < acceptable_rhodiff)
376 && (std::abs(env.p[i] / env.p[i + 1] - 1) < acceptable_pdiff)) {
384 const double rhomolar_vap_start = env.rhomolar_vap[i], rhomolar_vap_end = env.rhomolar_vap[i + 1];
386 double factor = pow(rhomolar_vap_end / rhomolar_vap_start, 1.0 / N);
388 int failure_count = 0;
389 for (
double rhomolar_vap = rhomolar_vap_start * factor; rhomolar_vap < rhomolar_vap_end; rhomolar_vap *= factor) {
391 IO.
x.resize(IO.
y.size());
392 if (i < env.T.size() - 3) {
395 for (std::size_t j = 0; j < IO.
x.size() - 1; ++j) {
401 for (std::size_t j = 0; j < IO.
x.size() - 1; ++j) {
405 IO.
x[IO.
x.size() - 1] = 1 - std::accumulate(IO.
x.begin(), IO.
x.end() - 1, 0.0);
407 NR.
call(HEOS, IO.
y, IO.
x, IO);
426 if (failure_count > 0) {
429 }
while (i < env.T.size() - 1);
432 int _i =
static_cast<int>(i);
433 std::vector<double>
const *x, *y;
443 y = &(env.rhomolar_vap);
446 y = &(env.hmolar_vap);
449 y = &(env.smolar_vap);
452 y = &(env.cpmolar_vap);
455 y = &(env.cvmolar_vap);
458 y = &(env.viscosity_vap);
461 y = &(env.conductivity_vap);
464 y = &(env.speed_sound_vap);
467 throw ValueError(
"Pointer to vector y is unset in is_inside");
470 double inval = value1;
480 x = &(env.rhomolar_vap);
483 x = &(env.hmolar_vap);
486 x = &(env.smolar_vap);
489 throw ValueError(
"Pointer to vector x is unset in is_inside");
491 if (_i + 2 >=
static_cast<int>(y->size())) {
494 if (_i + 1 >=
static_cast<int>(y->size())) {
501 double outval =
CubicInterp(*x, *y, _i - 1, _i, _i + 1, _i + 2, inval);
502 i =
static_cast<std::size_t
>(_i);
521 std::size_t iTmax = std::distance(env.T.begin(), std::max_element(env.T.begin(), env.T.end()));
524 std::size_t ipmax = std::distance(env.p.begin(), std::max_element(env.p.begin(), env.p.end()));
529 env.
TypeI = env.p[env.p.size() - 1] < env.p[ipmax];
534 for (
int imaxima = 0; imaxima <= 1; ++imaxima) {
535 maxima_points maxima;
536 if (imaxima == PMAX_SAT) {
538 }
else if (imaxima == TMAX_SAT) {
541 throw ValueError(
"I don't understand your maxima index");
546 if (maxima == TMAX_SAT) {
548 if (iTmax > env.T.size() - 3) {
551 spline.
add_4value_constraints(env.rhomolar_vap[iTmax - 1], env.rhomolar_vap[iTmax], env.rhomolar_vap[iTmax + 1],
552 env.rhomolar_vap[iTmax + 2], env.T[iTmax - 1], env.T[iTmax], env.T[iTmax + 1], env.T[iTmax + 2]);
555 if (ipmax > env.p.size() - 3) {
558 spline.
add_4value_constraints(env.rhomolar_vap[ipmax - 1], env.rhomolar_vap[ipmax], env.rhomolar_vap[ipmax + 1],
559 env.rhomolar_vap[ipmax + 2], env.p[ipmax - 1], env.p[ipmax], env.p[ipmax + 1], env.p[ipmax + 2]);
567 double rho0 = _HUGE, rho1 = _HUGE, rho2 = _HUGE;
568 solve_cubic(0, 3 * spline.
a, 2 * spline.
b, spline.
c, Nsoln, rho0, rho1, rho2);
575 }
else if (Nsoln == 2) {
576 if (
is_in_closed_range(env.rhomolar_vap[imax - 1], env.rhomolar_vap[imax + 1], rho0)) {
579 if (
is_in_closed_range(env.rhomolar_vap[imax - 1], env.rhomolar_vap[imax + 1], rho1)) {
583 throw ValueError(
"More than 2 solutions found");
590 maxima_points maxima;
595 this->HEOS = &HEOS, this->imax = imax;
596 this->maxima = maxima;
598 double call(
double rhomolar_vap) {
605 if (imax >= env.T.size() - 2) {
610 for (std::size_t i = 0; i < IO.
x.size() - 1; ++i)
614 IO.
x[IO.
x.size() - 1] = 1 - std::accumulate(IO.
x.begin(), IO.
x.end() - 1, 0.0);
615 NR.
call(*HEOS, IO.
y, IO.
x, IO);
616 if (maxima == TMAX_SAT) {
624 solver_resid resid(HEOS, imax, maxima);
630 if (rho > env.rhomolar_vap[imax]) {
634 env.
insert_variables(resid.IO.T, resid.IO.p, resid.IO.rhomolar_liq, resid.IO.rhomolar_vap, resid.IO.hmolar_liq, resid.IO.hmolar_vap,
635 resid.IO.smolar_liq, resid.IO.smolar_vap, resid.IO.x, resid.IO.y, imax);
643 env.
iTsat_max = std::distance(env.T.begin(), std::max_element(env.T.begin(), env.T.end()));
646 env.
ipsat_max = std::distance(env.p.begin(), std::max_element(env.p.begin(), env.p.end()));
651 std::vector<std::pair<std::size_t, std::size_t>> intersections;
653 for (std::size_t i = 0; i < env.p.size() - 1; ++i) {
654 bool matched =
false;
681 intersections.push_back(std::pair<std::size_t, std::size_t>(i, i + 1));
684 return intersections;
687 std::size_t& iclosest,
SimpleState& closest_state) {
689 std::vector<std::pair<std::size_t, std::size_t>> intersections =
find_intersections(env, iInput1, value1);
692 std::cout <<
format(
"is_inside(%Lg,%Lg); iTsat_max=%d; ipsat_max=%d\n", value1, value2, env.
iTsat_max, env.
ipsat_max);
703 if (intersections.size() == 0) {
704 throw ValueError(
format(
"Input is out of range for primary value [%Lg], inputs were (%s,%Lg,%s,%Lg); no intersections found", value1,
711 if (intersections.size() % 2 != 0) {
712 throw ValueError(
"Input is weird; odd number of intersections found");
716 if (intersections.size() % 2 == 0) {
717 if (intersections.size() != 2) {
718 throw ValueError(
"for now only even value accepted is 2");
720 std::vector<std::size_t> other_indices(4, 0);
721 std::vector<double>
const* y;
722 std::vector<double> other_values(4, 0);
723 other_indices[0] = intersections[0].first;
724 other_indices[1] = intersections[0].second;
725 other_indices[2] = intersections[1].first;
726 other_indices[3] = intersections[1].second;
736 y = &(env.rhomolar_vap);
739 y = &(env.hmolar_vap);
742 y = &(env.smolar_vap);
745 throw ValueError(
"Pointer to vector y is unset in is_inside");
748 other_values[0] = (*y)[other_indices[0]];
749 other_values[1] = (*y)[other_indices[1]];
750 other_values[2] = (*y)[other_indices[2]];
751 other_values[3] = (*y)[other_indices[3]];
753 CoolPropDbl min_other = *(std::min_element(other_values.begin(), other_values.end()));
754 CoolPropDbl max_other = *(std::max_element(other_values.begin(), other_values.end()));
757 std::cout <<
format(
"is_inside: min: %Lg max: %Lg val: %Lg\n", min_other, max_other, value2);
764 std::vector<CoolPropDbl> d(4, 0);
765 d[0] = std::abs(other_values[0] - value2);
766 d[1] = std::abs(other_values[1] - value2);
767 d[2] = std::abs(other_values[2] - value2);
768 d[3] = std::abs(other_values[3] - value2);
771 std::size_t idist = std::distance(d.begin(), std::min_element(d.begin(), d.end()));
773 iclosest = other_indices[idist];
778 closest_state.
T = env.T[iclosest];
779 closest_state.
p = env.p[iclosest];
780 closest_state.
rhomolar = env.rhomolar_vap[iclosest];
781 closest_state.
hmolar = env.hmolar_vap[iclosest];
782 closest_state.
smolar = env.smolar_vap[iclosest];
783 closest_state.
Q = env.Q[iclosest];
786 std::cout <<
format(
"is_inside: it is not inside") << std::endl;
796 if (std::abs(y1 - value2) < std::abs(y2 - value2)) {
797 iclosest = intersections[0].first;
799 iclosest = intersections[1].first;
804 closest_state.
T = env.T[iclosest];
805 closest_state.
p = env.p[iclosest];
806 closest_state.
rhomolar = env.rhomolar_vap[iclosest];
807 closest_state.
hmolar = env.hmolar_vap[iclosest];
808 closest_state.
smolar = env.smolar_vap[iclosest];
809 closest_state.
Q = env.Q[iclosest];
816 throw ValueError(
"You have a funny number of intersections in is_inside");