CoolProp  6.6.1dev
An open-source fluid property and humid air property database
PhaseEnvelopeRoutines.cpp
Go to the documentation of this file.
1 #ifndef PHASEENVELOPE_H
2 #define PHASEENVELOPE_H
3 
5 #include "VLERoutines.h"
7 #include "PhaseEnvelope.h"
8 #include "CoolPropTools.h"
9 #include "Configuration.h"
10 #include "CPnumerics.h"
11 
12 namespace CoolProp {
13 
14 void PhaseEnvelopeRoutines::build(HelmholtzEOSMixtureBackend& HEOS, const std::string& level) {
15  if (HEOS.get_mole_fractions_ref().empty()) {
16  throw ValueError("Mole fractions have not been set yet.");
17  }
18  bool debug = get_debug_level() > 0 || false;
19  if (HEOS.get_mole_fractions_ref().size() == 1) {
20  // It's a pure fluid
21  PhaseEnvelopeData& env = HEOS.PhaseEnvelope;
22  env.resize(HEOS.mole_fractions.size());
23 
24  // Breakpoints in the phase envelope
25  std::vector<CoolPropDbl> Tbp, Qbp;
26  std::vector<std::size_t> Nbp;
27  // Triple point vapor up to Tmax_sat
28  Tbp.push_back(HEOS.Ttriple());
29  Qbp.push_back(1.0);
30  Nbp.push_back(40);
31 
32  if (HEOS.is_pure()) {
33  // Up to critical point, back to triple point on the liquid side
34  Tbp.push_back(HEOS.T_critical() - 1e-3);
35  Qbp.push_back(0.0);
36  Tbp.push_back(HEOS.Ttriple());
37  Nbp.push_back(40);
38  } else {
39  SimpleState max_sat_T = HEOS.get_state("max_sat_T"), max_sat_p = HEOS.get_state("max_sat_p"), crit = HEOS.get_state("critical");
40  if (max_sat_T.rhomolar < crit.rhomolar && max_sat_T.rhomolar < max_sat_p.rhomolar) {
41  Tbp.push_back(HEOS.calc_Tmax_sat());
42  if (max_sat_p.rhomolar < crit.rhomolar) {
43  // psat_max density less than critical density
44  Qbp.push_back(1.0);
45  Qbp.push_back(1.0);
46  Tbp.push_back(max_sat_p.T);
47  Tbp.push_back(crit.T);
48  } else {
49  // Vapor line density less than critical density
50  Qbp.push_back(1.0);
51  Qbp.push_back(0.0);
52  Nbp.push_back(10);
53  Tbp.push_back(crit.T);
54  Tbp.push_back(max_sat_p.T);
55  }
56  Nbp.push_back(10);
57  Nbp.push_back(10);
58  Qbp.push_back(0.0);
59  Nbp.push_back(40);
60  Tbp.push_back(HEOS.Ttriple());
61  } else {
62  throw ValueError(format(""));
63  }
64  }
65 
66  for (std::size_t i = 0; i < Tbp.size() - 1; ++i) {
67  CoolPropDbl Tmin = Tbp[i], Tmax = Tbp[i + 1];
68  std::size_t N = Nbp[i];
69  for (CoolPropDbl T = Tmin; is_in_closed_range(Tmin, Tmax, T); T += (Tmax - Tmin) / (N - 1)) {
70  try {
71  HEOS.update(QT_INPUTS, Qbp[i], T);
72  } catch (...) {
73  continue;
74  }
75  if (Qbp[i] > 0.5) {
79  std::vector<CoolPropDbl>(1, 1.0), std::vector<CoolPropDbl>(1, 1.0));
80  } else {
84  std::vector<CoolPropDbl>(1, 1.0), std::vector<CoolPropDbl>(1, 1.0));
85  }
86  }
87  }
88  } else {
89  // It's a mixture
90  // --------------
91 
92  // First we try to generate all the critical points. This
93  // is very useful
94  std::vector<CriticalState> critpts;
95  // try{
96  // critpts = HEOS.all_critical_points();
97  // //throw CoolProp::ValueError("critical points disabled");
98  // }
99  // catch(std::exception &e)
100  // {
101  // if (debug){ std::cout << e.what() << std::endl; }
102  // };
103 
104  std::size_t failure_count = 0;
105  // Set some input options
108  io.Nstep_max = 20;
109 
110  // Set the pressure to a low pressure
111  HEOS._p = get_config_double(PHASE_ENVELOPE_STARTING_PRESSURE_PA); //[Pa]
112  HEOS._Q = 1;
113 
114  // Get an extremely rough guess by interpolation of ln(p) v. T curve where the limits are mole-fraction-weighted
116 
117  // Use Wilson iteration to obtain updated guess for temperature
119 
120  // Actually call the successive substitution solver
121  io.beta = 1;
122  SaturationSolvers::successive_substitution(HEOS, HEOS._Q, Tguess, HEOS._p, HEOS.mole_fractions, HEOS.K, io);
123 
124  // Use the residual function based on x_i, T and rho' as independent variables. rho'' is specified
127 
128  IO.bubble_point = false; // Do a "dewpoint" calculation all the way around
129  IO.x = io.x;
130  IO.y = HEOS.mole_fractions;
131  IO.rhomolar_liq = io.rhomolar_liq;
132  IO.rhomolar_vap = io.rhomolar_vap;
133  IO.T = io.T;
134  IO.p = io.p;
135  IO.Nstep_max = 30;
136 
137  /*
138  IO.p = 1e5;
139  IO.rhomolar_liq = 17257.17130;
140  IO.rhomolar_vap = 56.80022884;
141  IO.T = 219.5200523;
142  IO.x[0] = 0.6689704673;
143  IO.x[1] = 0.3310295327;
144  */
145 
146  //IO.rhomolar_liq *= 1.2;
147 
149 
150  NR.call(HEOS, IO.y, IO.x, IO);
151 
152  // Switch to density imposed
154 
155  bool dont_extrapolate = false;
156 
157  PhaseEnvelopeData& env = HEOS.PhaseEnvelope;
158  env.resize(HEOS.mole_fractions.size());
159 
160  std::size_t iter = 0, //< The iteration counter
161  iter0 = 0; //< A reference point for the counter, can be increased to go back to linear interpolation
162  CoolPropDbl factor = 1.05;
163 
164  for (;;) {
165  top_of_loop:; // A goto label so that nested loops can break out to the top of this loop
166 
167  if (failure_count > 5) {
168  // Stop since we are stuck at a bad point
169  //throw SolutionError("stuck");
170  return;
171  }
172 
173  if (iter - iter0 > 0) {
174  IO.rhomolar_vap *= factor;
175  }
176  if (dont_extrapolate) {
177  // Reset the step to a reasonably small size
178  factor = 1.0001;
179  } else if (iter - iter0 == 2) {
180  IO.T = LinearInterp(env.rhomolar_vap, env.T, iter - 2, iter - 1, IO.rhomolar_vap);
181  IO.rhomolar_liq = LinearInterp(env.rhomolar_vap, env.rhomolar_liq, iter - 2, iter - 1, IO.rhomolar_vap);
182  for (std::size_t i = 0; i < IO.x.size() - 1; ++i) // First N-1 elements
183  {
184  IO.x[i] = LinearInterp(env.rhomolar_vap, env.x[i], iter - 2, iter - 1, IO.rhomolar_vap);
185  }
186  } else if (iter - iter0 == 3) {
187  IO.T = QuadInterp(env.rhomolar_vap, env.T, iter - 3, iter - 2, iter - 1, IO.rhomolar_vap);
188  IO.rhomolar_liq = QuadInterp(env.rhomolar_vap, env.rhomolar_liq, iter - 3, iter - 2, iter - 1, IO.rhomolar_vap);
189  for (std::size_t i = 0; i < IO.x.size() - 1; ++i) // First N-1 elements
190  {
191  IO.x[i] = QuadInterp(env.rhomolar_vap, env.x[i], iter - 3, iter - 2, iter - 1, IO.rhomolar_vap);
192  }
193  } else if (iter - iter0 > 3) {
194  // Use the spline interpolation class of Devin Lane: http://shiftedbits.org/2011/01/30/cubic-spline-interpolation/
195  Spline<double, double> spl_T(env.rhomolar_vap, env.T);
196  IO.T = spl_T.interpolate(IO.rhomolar_vap);
197  Spline<double, double> spl_rho(env.rhomolar_vap, env.rhomolar_liq);
198  IO.rhomolar_liq = spl_rho.interpolate(IO.rhomolar_vap);
199 
200  // Check if there is a large deviation from linear interpolation - this suggests a step size that is so large that a minima or maxima of the interpolation function is crossed
201  CoolPropDbl T_linear = LinearInterp(env.rhomolar_vap, env.T, iter - 2, iter - 1, IO.rhomolar_vap);
202  if (std::abs((T_linear - IO.T) / IO.T) > 0.1) {
203  // Try again, but with a smaller step
204  IO.rhomolar_vap /= factor;
205  factor = 1 + (factor - 1) / 2;
206  failure_count++;
207  continue;
208  }
209  for (std::size_t i = 0; i < IO.x.size() - 1; ++i) // First N-1 elements
210  {
211  // Use the spline interpolation class of Devin Lane: http://shiftedbits.org/2011/01/30/cubic-spline-interpolation/
212  Spline<double, double> spl(env.rhomolar_vap, env.x[i]);
213  IO.x[i] = spl.interpolate(IO.rhomolar_vap);
214 
215  if (IO.x[i] < 0 || IO.x[i] > 1) {
216  // Try again, but with a smaller step
217  IO.rhomolar_vap /= factor;
218  factor = 1 + (factor - 1) / 2;
219  failure_count++;
220  goto top_of_loop;
221  }
222  }
223  }
224 
225  // The last mole fraction is sum of N-1 first elements
226  IO.x[IO.x.size() - 1] = 1 - std::accumulate(IO.x.begin(), IO.x.end() - 1, 0.0);
227 
228  // Uncomment to check guess values for Newton-Raphson
229  //std::cout << "\t\tdv " << IO.rhomolar_vap << " dl " << IO.rhomolar_liq << " T " << IO.T << " x " << vec_to_string(IO.x, "%0.10Lg") << std::endl;
230 
231  // Dewpoint calculation, liquid (x) is incipient phase
232  try {
233  NR.call(HEOS, IO.y, IO.x, IO);
234  if (!ValidNumber(IO.rhomolar_liq) || !ValidNumber(IO.p) || !ValidNumber(IO.T)) {
235  throw ValueError("Invalid number");
236  }
237  // Reject trivial solution
238  if (std::abs(IO.rhomolar_liq - IO.rhomolar_vap) < 1e-3) {
239  throw ValueError("Trivial solution");
240  }
241  // Reject negative presssure
242  if (IO.p < 0) {
243  throw ValueError("negative pressure");
244  }
245  // Reject steps with enormous steps in temperature
246  if (!env.T.empty() && std::abs(env.T[env.T.size() - 1] - IO.T) > 100) {
247  throw ValueError("Change in temperature too large");
248  }
249  } catch (std::exception& e) {
250  if (debug) {
251  std::cout << e.what() << std::endl;
252  }
253  //std::cout << IO.T << " " << IO.p << std::endl;
254  // Try again, but with a smaller step
255  IO.rhomolar_vap /= factor;
256  if (iter < 4) {
257  throw ValueError(format("Unable to calculate at least 4 points in phase envelope; quitting"));
258  }
259  IO.rhomolar_liq = QuadInterp(env.rhomolar_vap, env.rhomolar_liq, iter - 3, iter - 2, iter - 1, IO.rhomolar_vap);
260  factor = 1 + (factor - 1) / 2;
261  failure_count++;
262  continue;
263  }
264 
265  if (debug) {
266  std::cout << "dv " << IO.rhomolar_vap << " dl " << IO.rhomolar_liq << " T " << IO.T << " p " << IO.p << " hl " << IO.hmolar_liq
267  << " hv " << IO.hmolar_vap << " sl " << IO.smolar_liq << " sv " << IO.smolar_vap << " x " << vec_to_string(IO.x, "%0.10Lg")
268  << " Ns " << IO.Nsteps << " factor " << factor << std::endl;
269  }
270  env.store_variables(IO.T, IO.p, IO.rhomolar_liq, IO.rhomolar_vap, IO.hmolar_liq, IO.hmolar_vap, IO.smolar_liq, IO.smolar_vap, IO.x, IO.y);
271 
272  iter++;
273 
274  // CoolPropDbl abs_rho_difference = std::abs((IO.rhomolar_liq - IO.rhomolar_vap)/IO.rhomolar_liq);
275 
276  // bool next_crosses_crit = false;
277  // if (it_critpts != critpts.end() ){
278  // // Density at the next critical point
279  // double rhoc = (*it_critpts).rhomolar;
280  // // Next vapor density that will be used
281  // double rho_next = IO.rhomolar_vap*factor;
282  // // If the signs of the differences are different, you have crossed
283  // // the critical point density and have a phase inversion
284  // // on your hands
285  // next_crosses_crit = ((IO.rhomolar_vap-rhoc)*(rho_next-rhoc) < 0);
286  // }
287 
288  // // Critical point jump
289  // if (next_crosses_crit || (abs_rho_difference < 0.01 && IO.rhomolar_liq > IO.rhomolar_vap)){
290  // //std::cout << "dv" << IO.rhomolar_vap << " dl " << IO.rhomolar_liq << " " << vec_to_string(IO.x, "%0.10Lg") << " " << vec_to_string(IO.y, "%0.10Lg") << std::endl;
291  // CoolPropDbl rhoc_approx = 0.5*IO.rhomolar_liq + 0.5*IO.rhomolar_vap;
292  // if (it_critpts != critpts.end() ){
293  // // We actually know what the critical point is to numerical precision
294  // rhoc_approx = (*it_critpts).rhomolar;
295  // }
296  // CoolPropDbl rho_vap_new = 1.05*rhoc_approx;
297  // // Linearly interpolate to get new guess for T
298  // IO.T = LinearInterp(env.rhomolar_vap,env.T,iter-2,iter-1,rho_vap_new);
299  // IO.rhomolar_liq = LinearInterp(env.rhomolar_vap, env.rhomolar_liq, iter-2, iter-1, rho_vap_new);
300  // for (std::size_t i = 0; i < IO.x.size()-1; ++i){
301  // IO.x[i] = CubicInterp(env.rhomolar_vap, env.x[i], iter-4, iter-3, iter-2, iter-1, rho_vap_new);
302  // }
303  // IO.x[IO.x.size()-1] = 1 - std::accumulate(IO.x.begin(), IO.x.end()-1, 0.0);
304  // factor = rho_vap_new/IO.rhomolar_vap;
305  // dont_extrapolate = true; // So that we use the mole fractions we calculated here instead of the extrapolated values
306  // if (debug) std::cout << "[CRIT jump] new values: dv " << rho_vap_new << " dl " << IO.rhomolar_liq << " " << vec_to_string(IO.x, "%0.10Lg") << " " << vec_to_string(IO.y, "%0.10Lg") << std::endl;
307  // iter0 = iter - 1; // Back to linear interpolation again
308  // continue;
309  // }
310 
311  dont_extrapolate = false;
312  if (iter < 5) {
313  continue;
314  }
315  if (IO.Nsteps > 10) {
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;
321  }
322  // Min step is 1.01
323  factor = std::max(factor, static_cast<CoolPropDbl>(1.01));
324  // As we approach the critical point, control step size
325  if (std::abs(IO.rhomolar_liq / IO.rhomolar_vap - 1) < 4) {
326  // Max step is 1.1
327  factor = std::min(factor, static_cast<CoolPropDbl>(1.1));
328  }
329 
330  // Stop if the pressure is below the starting pressure
331  // or if the composition of one of the phases becomes almost pure
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)) {
334  env.built = true;
335  if (debug) {
336  std::cout << format("envelope built.\n");
337  std::cout << format("closest fraction to 1.0: distance %g\n", 1 - max_fraction);
338  }
339 
340  // Now we refine the phase envelope to add some points in places that are still pretty rough
341  refine(HEOS, level);
342 
343  return;
344  }
345 
346  // Reset the failure counter
347  failure_count = 0;
348  }
349  }
350 }
351 
352 void PhaseEnvelopeRoutines::refine(HelmholtzEOSMixtureBackend& HEOS, const std::string& level) {
353  bool debug = (get_debug_level() > 0 || false);
354  PhaseEnvelopeData& env = HEOS.PhaseEnvelope;
358  IO.bubble_point = false;
359  IO.y = HEOS.get_mole_fractions();
360 
361  double acceptable_pdiff = 0.5;
362  double acceptable_rhodiff = 0.25;
363  int N = 5; // Number of steps of refining
364  if (level == "veryfine") {
365  acceptable_pdiff = 0.1;
366  acceptable_rhodiff = 0.1;
367  }
368  if (level == "none") {
369  return;
370  }
371  std::size_t i = 0;
372  do {
373 
374  // Don't do anything if change in density and pressure is small enough
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)) {
377  i++;
378  continue;
379  }
380 
381  // Ok, now we are going to do some more refining in this step
382 
383  // Vapor densities for this step, vapor density monotonically increasing
384  const double rhomolar_vap_start = env.rhomolar_vap[i], rhomolar_vap_end = env.rhomolar_vap[i + 1];
385 
386  double factor = pow(rhomolar_vap_end / rhomolar_vap_start, 1.0 / N);
387 
388  int failure_count = 0;
389  for (double rhomolar_vap = rhomolar_vap_start * factor; rhomolar_vap < rhomolar_vap_end; rhomolar_vap *= factor) {
390  IO.rhomolar_vap = rhomolar_vap;
391  IO.x.resize(IO.y.size());
392  if (i < env.T.size() - 3) {
393  IO.T = CubicInterp(env.rhomolar_vap, env.T, i, i + 1, i + 2, i + 3, IO.rhomolar_vap);
394  IO.rhomolar_liq = CubicInterp(env.rhomolar_vap, env.rhomolar_liq, i, i + 1, i + 2, i + 3, IO.rhomolar_vap);
395  for (std::size_t j = 0; j < IO.x.size() - 1; ++j) { // First N-1 elements
396  IO.x[j] = CubicInterp(env.rhomolar_vap, env.x[j], i, i + 1, i + 2, i + 3, IO.rhomolar_vap);
397  }
398  } else {
399  IO.T = CubicInterp(env.rhomolar_vap, env.T, i, i - 1, i - 2, i - 3, IO.rhomolar_vap);
400  IO.rhomolar_liq = CubicInterp(env.rhomolar_vap, env.rhomolar_liq, i, i - 1, i - 2, i - 3, IO.rhomolar_vap);
401  for (std::size_t j = 0; j < IO.x.size() - 1; ++j) { // First N-1 elements
402  IO.x[j] = CubicInterp(env.rhomolar_vap, env.x[j], i, i - 1, i - 2, i - 3, IO.rhomolar_vap);
403  }
404  }
405  IO.x[IO.x.size() - 1] = 1 - std::accumulate(IO.x.begin(), IO.x.end() - 1, 0.0);
406  try {
407  NR.call(HEOS, IO.y, IO.x, IO);
408  if (!ValidNumber(IO.rhomolar_liq) || !ValidNumber(IO.p)) {
409  throw ValueError("invalid numbers");
410  }
411  env.insert_variables(IO.T, IO.p, IO.rhomolar_liq, IO.rhomolar_vap, IO.hmolar_liq, IO.hmolar_vap, IO.smolar_liq, IO.smolar_vap, IO.x,
412  IO.y, i + 1);
413  if (debug) {
414  std::cout << "dv " << IO.rhomolar_vap << " dl " << IO.rhomolar_liq << " T " << IO.T << " p " << IO.p << " hl " << IO.hmolar_liq
415  << " hv " << IO.hmolar_vap << " sl " << IO.smolar_liq << " sv " << IO.smolar_vap << " x "
416  << vec_to_string(IO.x, "%0.10Lg") << " Ns " << IO.Nsteps << std::endl;
417  }
418  } catch (...) {
419  failure_count++;
420  continue;
421  }
422  i++;
423  }
424  // If we had a failure, we don't want to get stuck on this value of i,
425  // so we bump up one and keep moving
426  if (failure_count > 0) {
427  i++;
428  }
429  } while (i < env.T.size() - 1);
430 }
431 double PhaseEnvelopeRoutines::evaluate(const PhaseEnvelopeData& env, parameters output, parameters iInput1, double value1, std::size_t& i) {
432  int _i = static_cast<int>(i);
433  std::vector<double> const *x, *y;
434 
435  switch (output) {
436  case iT:
437  y = &(env.T);
438  break;
439  case iP:
440  y = &(env.p);
441  break;
442  case iDmolar:
443  y = &(env.rhomolar_vap);
444  break;
445  case iHmolar:
446  y = &(env.hmolar_vap);
447  break;
448  case iSmolar:
449  y = &(env.smolar_vap);
450  break;
451  case iCpmolar:
452  y = &(env.cpmolar_vap);
453  break;
454  case iCvmolar:
455  y = &(env.cvmolar_vap);
456  break;
457  case iviscosity:
458  y = &(env.viscosity_vap);
459  break;
460  case iconductivity:
461  y = &(env.conductivity_vap);
462  break;
463  case ispeed_sound:
464  y = &(env.speed_sound_vap);
465  break;
466  default:
467  throw ValueError("Pointer to vector y is unset in is_inside");
468  }
469 
470  double inval = value1;
471  switch (iInput1) {
472  case iT:
473  x = &(env.T);
474  break;
475  case iP:
476  x = &(env.lnp);
477  inval = log(value1);
478  break;
479  case iDmolar:
480  x = &(env.rhomolar_vap);
481  break;
482  case iHmolar:
483  x = &(env.hmolar_vap);
484  break;
485  case iSmolar:
486  x = &(env.smolar_vap);
487  break;
488  default:
489  throw ValueError("Pointer to vector x is unset in is_inside");
490  }
491  if (_i + 2 >= static_cast<int>(y->size())) {
492  _i--;
493  }
494  if (_i + 1 >= static_cast<int>(y->size())) {
495  _i--;
496  }
497  if (_i - 1 < 0) {
498  _i++;
499  }
500 
501  double outval = CubicInterp(*x, *y, _i - 1, _i, _i + 1, _i + 2, inval);
502  i = static_cast<std::size_t>(_i);
503  return outval;
504 }
506  // No finalization for pure or pseudo-pure fluids
507  if (HEOS.get_mole_fractions_ref().size() == 1) {
508  return;
509  }
510 
511  enum maxima_points
512  {
513  PMAX_SAT = 0,
514  TMAX_SAT = 1
515  };
516  std::size_t imax; // Index of the maximal temperature or pressure
517 
518  PhaseEnvelopeData& env = HEOS.PhaseEnvelope;
519 
520  // Find the index of the point with the highest temperature
521  std::size_t iTmax = std::distance(env.T.begin(), std::max_element(env.T.begin(), env.T.end()));
522 
523  // Find the index of the point with the highest pressure
524  std::size_t ipmax = std::distance(env.p.begin(), std::max_element(env.p.begin(), env.p.end()));
525 
526  // Determine if the phase envelope corresponds to a Type I mixture
527  // For now we consider a mixture to be Type I if the pressure at the
528  // end of the envelope is lower than max pressure pressure
529  env.TypeI = env.p[env.p.size() - 1] < env.p[ipmax];
530 
531  // Approximate solutions for the maxima of the phase envelope
532  // See method in Gernert. We use our spline class to find the coefficients
533  if (env.TypeI) {
534  for (int imaxima = 0; imaxima <= 1; ++imaxima) {
535  maxima_points maxima;
536  if (imaxima == PMAX_SAT) {
537  maxima = PMAX_SAT;
538  } else if (imaxima == TMAX_SAT) {
539  maxima = TMAX_SAT;
540  } else {
541  throw ValueError("I don't understand your maxima index");
542  }
543 
544  // Spline using the points around it
545  SplineClass spline;
546  if (maxima == TMAX_SAT) {
547  imax = iTmax;
548  if (iTmax > env.T.size() - 3) {
549  iTmax -= 2;
550  }
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]);
553  } else {
554  imax = ipmax;
555  if (ipmax > env.p.size() - 3) {
556  ipmax -= 2;
557  }
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]);
560  }
561  spline.build(); // y = a*rho^3 + b*rho^2 + c*rho + d
562 
563  // Take derivative
564  // dy/drho = 3*a*rho^2 + 2*b*rho + c
565  // Solve quadratic for derivative to find rho
566  int Nsoln = 0;
567  double rho0 = _HUGE, rho1 = _HUGE, rho2 = _HUGE;
568  solve_cubic(0, 3 * spline.a, 2 * spline.b, spline.c, Nsoln, rho0, rho1, rho2);
569 
571  IO.rhomolar_vap = _HUGE;
572  // Find the correct solution
573  if (Nsoln == 1) {
574  IO.rhomolar_vap = rho0;
575  } else if (Nsoln == 2) {
576  if (is_in_closed_range(env.rhomolar_vap[imax - 1], env.rhomolar_vap[imax + 1], rho0)) {
577  IO.rhomolar_vap = rho0;
578  }
579  if (is_in_closed_range(env.rhomolar_vap[imax - 1], env.rhomolar_vap[imax + 1], rho1)) {
580  IO.rhomolar_vap = rho1;
581  }
582  } else {
583  throw ValueError("More than 2 solutions found");
584  }
585 
586  class solver_resid : public FuncWrapper1D
587  {
588  public:
589  std::size_t imax;
590  maxima_points maxima;
594  solver_resid(HelmholtzEOSMixtureBackend& HEOS, std::size_t imax, maxima_points maxima) {
595  this->HEOS = &HEOS, this->imax = imax;
596  this->maxima = maxima;
597  };
598  double call(double rhomolar_vap) {
599  PhaseEnvelopeData& env = HEOS->PhaseEnvelope;
601  IO.bubble_point = false;
602  IO.rhomolar_vap = rhomolar_vap;
603  IO.y = HEOS->get_mole_fractions();
604  IO.x = IO.y; // Just to give it good size
605  if (imax >= env.T.size() - 2) {
606  imax -= 2;
607  }
608  IO.T = CubicInterp(env.rhomolar_vap, env.T, imax - 1, imax, imax + 1, imax + 2, IO.rhomolar_vap);
609  IO.rhomolar_liq = CubicInterp(env.rhomolar_vap, env.rhomolar_liq, imax - 1, imax, imax + 1, imax + 2, IO.rhomolar_vap);
610  for (std::size_t i = 0; i < IO.x.size() - 1; ++i) // First N-1 elements
611  {
612  IO.x[i] = CubicInterp(env.rhomolar_vap, env.x[i], imax - 1, imax, imax + 1, imax + 2, IO.rhomolar_vap);
613  }
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) {
617  return NR.dTsat_dPsat;
618  } else {
619  return NR.dPsat_dTsat;
620  }
621  };
622  };
623 
624  solver_resid resid(HEOS, imax, maxima);
625  try {
626  double rho = Brent(resid, IO.rhomolar_vap * 0.95, IO.rhomolar_vap * 1.05, DBL_EPSILON, 1e-12, 100);
627 
628  // If maxima point is greater than density at point from the phase envelope, increase index by 1 so that the
629  // insertion will happen *after* the point in the envelope since density is monotonically increasing.
630  if (rho > env.rhomolar_vap[imax]) {
631  imax++;
632  }
633 
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);
636  } catch (...) {
637  // Don't do the insertion
638  }
639  }
640  }
641 
642  // Find the index of the point with the highest temperature
643  env.iTsat_max = std::distance(env.T.begin(), std::max_element(env.T.begin(), env.T.end()));
644 
645  // Find the index of the point with the highest pressure
646  env.ipsat_max = std::distance(env.p.begin(), std::max_element(env.p.begin(), env.p.end()));
647 }
648 
649 std::vector<std::pair<std::size_t, std::size_t>> PhaseEnvelopeRoutines::find_intersections(const PhaseEnvelopeData& env, parameters iInput,
650  double value) {
651  std::vector<std::pair<std::size_t, std::size_t>> intersections;
652 
653  for (std::size_t i = 0; i < env.p.size() - 1; ++i) {
654  bool matched = false;
655  switch (iInput) {
656  case iP:
657  if (is_in_closed_range(env.p[i], env.p[i + 1], value)) {
658  matched = true;
659  }
660  break;
661  case iT:
662  if (is_in_closed_range(env.T[i], env.T[i + 1], value)) {
663  matched = true;
664  }
665  break;
666  case iHmolar:
667  if (is_in_closed_range(env.hmolar_vap[i], env.hmolar_vap[i + 1], value)) {
668  matched = true;
669  }
670  break;
671  case iSmolar:
672  if (is_in_closed_range(env.smolar_vap[i], env.smolar_vap[i + 1], value)) {
673  matched = true;
674  }
675  break;
676  default:
677  throw ValueError(format("bad index to find_intersections"));
678  }
679 
680  if (matched) {
681  intersections.push_back(std::pair<std::size_t, std::size_t>(i, i + 1));
682  }
683  }
684  return intersections;
685 }
687  std::size_t& iclosest, SimpleState& closest_state) {
688  // Find the indices that bound the solution(s)
689  std::vector<std::pair<std::size_t, std::size_t>> intersections = find_intersections(env, iInput1, value1);
690 
691  if (get_debug_level() > 5) {
692  std::cout << format("is_inside(%Lg,%Lg); iTsat_max=%d; ipsat_max=%d\n", value1, value2, env.iTsat_max, env.ipsat_max);
693  }
694  // Check whether input is above max value
695  if (iInput1 == iT && 0 < env.iTsat_max && env.iTsat_max < env.T.size() && value1 > env.T[env.iTsat_max]) {
696  return false;
697  }
698  if (iInput1 == iP && 0 < env.ipsat_max && env.ipsat_max < env.p.size() && value1 > env.p[env.ipsat_max]) {
699  return false;
700  }
701 
702  // If number of intersections is 0, input is out of range, quit
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,
705  get_parameter_information(iInput1, "short").c_str(), value1, get_parameter_information(iInput2, "short").c_str(),
706  value2));
707  }
708 
709  // If number of intersections is 1, input will be determined based on the single intersection
710  // Need to know if values increase or decrease to the right of the intersection point
711  if (intersections.size() % 2 != 0) {
712  throw ValueError("Input is weird; odd number of intersections found");
713  }
714 
715  // If number of intersections is even, might be a bound
716  if (intersections.size() % 2 == 0) {
717  if (intersections.size() != 2) {
718  throw ValueError("for now only even value accepted is 2");
719  }
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;
727 
728  switch (iInput2) {
729  case iT:
730  y = &(env.T);
731  break;
732  case iP:
733  y = &(env.p);
734  break;
735  case iDmolar:
736  y = &(env.rhomolar_vap);
737  break;
738  case iHmolar:
739  y = &(env.hmolar_vap);
740  break;
741  case iSmolar:
742  y = &(env.smolar_vap);
743  break;
744  default:
745  throw ValueError("Pointer to vector y is unset in is_inside");
746  }
747 
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]];
752 
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()));
755 
756  if (get_debug_level() > 5) {
757  std::cout << format("is_inside: min: %Lg max: %Lg val: %Lg\n", min_other, max_other, value2);
758  }
759 
760  // If by using the outer bounds of the second variable, we are outside the range,
761  // then the value is definitely not inside the phase envelope and we don't need to
762  // do any more analysis.
763  if (!is_in_closed_range(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);
769 
770  // Index of minimum distance in the other_values vector
771  std::size_t idist = std::distance(d.begin(), std::min_element(d.begin(), d.end()));
772  // Index of closest point in the phase envelope
773  iclosest = other_indices[idist];
774 
775  // Get the state for the point which is closest to the desired value - this
776  // can be used as a bounding value in the outer single-phase flash routine
777  // since you know (100%) that it is a good bound
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];
784 
785  if (get_debug_level() > 5) {
786  std::cout << format("is_inside: it is not inside") << std::endl;
787  }
788  return false;
789  } else {
790  // Now we have to do a saturation flash call in order to determine whether or not we are inside the phase envelope or not
791 
792  // First we can interpolate using the phase envelope to get good guesses for the necessary values
793  CoolPropDbl y1 = evaluate(env, iInput2, iInput1, value1, intersections[0].first);
794  CoolPropDbl y2 = evaluate(env, iInput2, iInput1, value1, intersections[1].first);
795  if (is_in_closed_range(y1, y2, value2)) {
796  if (std::abs(y1 - value2) < std::abs(y2 - value2)) {
797  iclosest = intersections[0].first;
798  } else {
799  iclosest = intersections[1].first;
800  }
801  // Get the state for the point which is closest to the desired value - this
802  // can be used as a bounding value in the outer single-phase flash routine
803  // since you know (100%) that it is a good bound
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];
810  return true;
811  } else {
812  return false;
813  }
814  }
815  } else {
816  throw ValueError("You have a funny number of intersections in is_inside");
817  }
818 }
819 
820 } /* namespace CoolProp */
821 
822 #endif