Loading [MathJax]/extensions/TeX/AMSsymbols.js
CoolProp  6.7.1dev
An open-source fluid property and humid air property database
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
CoolPropPlot.cpp
Go to the documentation of this file.
1 #include "CoolPropPlot.h"
2 #include "CoolProp.h"
3 #include "CPnumerics.h"
4 #include <map>
5 
6 namespace CoolProp {
7 namespace Plot {
8 
9 namespace Detail {
10 
11 const double NaN = std::numeric_limits<double>::quiet_NaN();
12 
13 const int TS = CoolProp::iT * 10 + CoolProp::iSmass;
14 const int PH = CoolProp::iP * 10 + CoolProp::iHmass;
15 const int HS = CoolProp::iHmass * 10 + CoolProp::iSmass;
16 const int PS = CoolProp::iP * 10 + CoolProp::iSmass;
17 const int PD = CoolProp::iP * 10 + CoolProp::iDmass;
18 const int TD = CoolProp::iT * 10 + CoolProp::iDmass;
19 const int PT = CoolProp::iP * 10 + CoolProp::iT;
20 
21 enum IsolineSupported
22 {
23  No = 0,
24  Yes = 1,
25  Flipped = 2
26 };
27 
28 const std::map<CoolProp::parameters, std::map<int, IsolineSupported>> xy_switch = {
29  {CoolProp::iDmass, {{TS, Flipped}, {PH, Flipped}, {HS, Yes }, {PS, Flipped}, {PD, No }, {TD, No }, {PT, Yes }}},
30  {CoolProp::iHmass, {{TS, Yes }, {PH, No }, {HS, No }, {PS, Flipped}, {PD, Flipped}, {TD, Yes }, {PT, Yes }}},
31  {CoolProp::iP, {{TS, Yes }, {PH, No }, {HS, Yes }, {PS, No }, {PD, No }, {TD, Yes }, {PT, No }}},
32  {CoolProp::iSmass, {{TS, No }, {PH, Flipped}, {HS, No }, {PS, No }, {PD, Flipped}, {TD, Yes }, {PT, Flipped}}},
33  {CoolProp::iT, {{TS, No }, {PH, Flipped}, {HS, Yes }, {PS, Yes }, {PD, Yes }, {TD, No }, {PT, No }}},
34  {CoolProp::iQ, {{TS, Flipped}, {PH, Flipped}, {HS, Flipped}, {PS, Flipped}, {PD, Flipped}, {TD, Flipped}, {PT, Yes }}}
35 };
36 
37 Scale default_scale(CoolProp::parameters key) {
38  switch (key) {
39  case CoolProp::iDmass: return Scale::Log;
40  case CoolProp::iHmass: return Scale::Lin;
41  case CoolProp::iP: return Scale::Log;
42  case CoolProp::iSmass: return Scale::Lin;
43  case CoolProp::iT: return Scale::Lin;
44  case CoolProp::iUmass: return Scale::Lin;
45  case CoolProp::iQ: return Scale::Lin;
46  default: return Scale::Lin;
47  }
48 }
49 
50 inline std::shared_ptr<CoolProp::AbstractState> process_fluid_state(const std::string& fluid_ref) {
51  std::string backend;
52  std::string fluids;
53  CoolProp::extract_backend(fluid_ref, backend, fluids);
54  std::vector<double> fractions;
55  fluids = CoolProp::extract_fractions(fluids, fractions);
56 
57  return std::shared_ptr<CoolProp::AbstractState>(CoolProp::AbstractState::factory(backend, fluids));
58 }
59 
60 std::shared_ptr<CoolProp::AbstractState> get_critical_point(const std::shared_ptr<CoolProp::AbstractState>& state) {
61  CoolProp::CriticalState crit_state;
62  crit_state.T = Detail::NaN;
63  crit_state.p = Detail::NaN;
64  crit_state.rhomolar = Detail::NaN;
65  crit_state.rhomolar = Detail::NaN;
66  crit_state.stable = false;
67  try {
68  crit_state.T = state->T_critical();
69  crit_state.p = state->p_critical();
70  crit_state.rhomolar = state->rhomolar_critical();
71  crit_state.stable = true;
72  } catch (...) {
73  try {
74  for (CoolProp::CriticalState crit_state_tmp: state->all_critical_points()) {
75  if (crit_state_tmp.stable && (crit_state_tmp.T > crit_state.T || !std::isfinite(crit_state.T))) {
76  crit_state.T = crit_state_tmp.T;
77  crit_state.p = crit_state_tmp.p;
78  crit_state.rhomolar = crit_state_tmp.rhomolar;
79  crit_state.stable = crit_state_tmp.stable;
80  }
81  }
82  } catch (...) {
83  throw CoolProp::ValueError("Could not calculate the critical point data.");
84  }
85  }
86 
87  std::shared_ptr<CoolProp::AbstractState> new_state(CoolProp::AbstractState::factory(state->backend_name(), state->fluid_names()));
88  std::vector<double> masses = state->get_mass_fractions();
89  if (masses.size() > 1)
90  new_state->set_mass_fractions(masses);
91 
92  if (std::isfinite(crit_state.p) && std::isfinite(crit_state.T)) {
93  try {
94  new_state->specify_phase(CoolProp::iphase_critical_point);
95  new_state->update(CoolProp::PT_INPUTS, crit_state.p, crit_state.T);
96  return new_state;
97  } catch (...) { }
98  try {
99  new_state->update(CoolProp::PT_INPUTS, crit_state.p, crit_state.T);
100  return new_state;
101  } catch (...) { }
102  }
103 
104  if (std::isfinite(crit_state.rhomolar) && std::isfinite(crit_state.T)) {
105  try {
106  new_state->specify_phase(CoolProp::iphase_critical_point);
107  new_state->update(CoolProp::DmolarT_INPUTS, crit_state.rhomolar, crit_state.T);
108  return new_state;
109  } catch (...) { }
110  try {
111  new_state->update(CoolProp::DmolarT_INPUTS, crit_state.rhomolar, crit_state.T);
112  return new_state;
113  } catch (...) { }
114  }
115  throw CoolProp::ValueError("Could not calculate the critical point data.");
116  return nullptr;
117 }
118 
119 } /* namespace Detail */
120 
121 
122 std::vector<double> generate_values_in_range(Scale scale, const Range& range, int count) {
123  if (scale == Scale::Log)
124  return logspace(range.min, range.max, count);
125  else
126  return linspace(range.min, range.max, count);
127 }
128 
129 std::vector<double> generate_values_in_range(CoolProp::parameters type, const Range& range, int count) {
130  return generate_values_in_range(Detail::default_scale(type), range, count);
131 }
132 
133 
134 Isoline::Isoline(CoolProp::parameters key, CoolProp::parameters xkey, CoolProp::parameters ykey, double value, const std::shared_ptr<CoolProp::AbstractState>& state)
135  : key_(key),
136  xkey_(xkey),
137  ykey_(ykey),
138  value(value),
139  state_(state) {
140  this->critical_state_ = Detail::get_critical_point(state);
141 }
142 
143 Range Isoline::get_sat_bounds(CoolProp::parameters key) const {
144  double s = 1e-7;
145  double t_small = critical_state_->keyed_output(CoolProp::iT) * s;
146  double p_small = critical_state_->keyed_output(CoolProp::iP) * s;
147 
148  double t_triple = state_->trivial_keyed_output(CoolProp::iT_triple);
149  double t_min;
150  try {
151  t_min = state_->trivial_keyed_output(CoolProp::iT_min);
152  } catch (...) {
153  t_min = t_triple;
154  }
155  state_->update(CoolProp::QT_INPUTS, 0, std::max(t_triple, t_min) + t_small);
156  if (key == CoolProp::iP)
157  return {state_->keyed_output(CoolProp::iP) + p_small, critical_state_->keyed_output(CoolProp::iP) - p_small};
158  else if (key == CoolProp::iT)
159  return {state_->keyed_output(CoolProp::iT) + t_small, critical_state_->keyed_output(CoolProp::iT) - t_small};
160  else
161  throw CoolProp::ValueError("Invalid key");
162 }
163 
164 void Isoline::calc_sat_range(int count) {
165  Range t = get_sat_bounds(CoolProp::iT);
166  std::vector<double> two = ::linspace(t.min, t.max, count);
167  std::vector<double> one(two.size(), value);
168  CoolProp::input_pairs input_pair = CoolProp::QT_INPUTS;
169 
170  double t_crit = critical_state_->keyed_output(CoolProp::iT);
171  double p_crit = critical_state_->keyed_output(CoolProp::iP);
172  double x_crit = critical_state_->keyed_output(xkey_);
173  double y_crit = critical_state_->keyed_output(ykey_);
174  x.resize(one.size());
175  y.resize(one.size());
176  for (int i = 0; i < one.size(); ++i) {
177  try {
178  state_->update(input_pair, one[i], two[i]);
179  x[i] = state_->keyed_output(xkey_);
180  y[i] = state_->keyed_output(ykey_);
181  } catch (...) {
182  if (input_pair == CoolProp::QT_INPUTS && abs(two[i] - t_crit) < 1e0
183  || input_pair == CoolProp::PQ_INPUTS && abs(one[i] - p_crit) < 1e2) {
184  x[i] = x_crit;
185  y[i] = y_crit;
186  std::cerr << "ERROR near critical inputs" << std::endl;
187  } else {
188  x[i] = Detail::NaN;
189  y[i] = Detail::NaN;
190  std::cerr << "ERROR" << std::endl;
191  }
192  }
193  }
194 }
195 
196 void Isoline::update_pair(int& ipos, int& xpos, int& ypos, int& pair) {
197  Detail::IsolineSupported should_switch = Detail::xy_switch.at(key_).at(ykey_ * 10 + xkey_);
198  double out1, out2;
199  if (should_switch == Detail::IsolineSupported::No)
200  throw CoolProp::ValueError("This isoline cannot be calculated!");
201  else if (should_switch == Detail::IsolineSupported::Yes)
202  pair = CoolProp::generate_update_pair(key_, 0.0, xkey_, 1.0, out1, out2);
203  else if (should_switch == Detail::IsolineSupported::Flipped)
204  pair = CoolProp::generate_update_pair(key_, 0.0, ykey_, 1.0, out1, out2);
205 
206  bool should_swap = (out1 != 0.0);
207 
208  if (should_switch == Detail::IsolineSupported::Yes && !should_swap) {
209  ipos = 0;
210  xpos = 1;
211  ypos = 2;
212  } else if (should_switch == Detail::IsolineSupported::Flipped && !should_swap) {
213  ipos = 0;
214  xpos = 2;
215  ypos = 1;
216  } else if (should_switch == Detail::IsolineSupported::Yes && should_swap) {
217  ipos = 1;
218  xpos = 0;
219  ypos = 2;
220  } else if (should_switch == Detail::IsolineSupported::Flipped && should_swap) {
221  ipos = 1;
222  xpos = 2;
223  ypos = 0;
224  } else {
225  throw CoolProp::ValueError("Check the code, this should not happen!");
226  }
227 }
228 
229 void Isoline::calc_range(std::vector<double>& xvals, std::vector<double>& yvals) {
230  if (key_ == CoolProp::iQ) {
231  calc_sat_range(xvals.size());
232  } else {
233  int ipos, xpos, ypos, pair;
234  update_pair(ipos, xpos, ypos, pair);
235 
236  std::vector<double> ivals(xvals.size(), value);
237  std::vector<int> order = {ipos, xpos, ypos};
238  std::vector<CoolProp::parameters> idxs(3);
239  idxs[ipos] = key_;
240  idxs[xpos] = xkey_;
241  idxs[ypos] = ykey_;
242  std::vector<std::vector<double>> vals(3);
243  vals[ipos] = ivals;
244  vals[xpos] = xvals;
245  vals[ypos] = yvals;
246 
247  for (int i = 0; i < vals[2].size(); ++i) {
248  try {
249  state_->update((CoolProp::input_pairs)pair, vals[0][i], vals[1][i]);
250  vals[2][i] = state_->keyed_output(idxs[2]);
251  } catch (...) {
252  vals[2][i] = Detail::NaN;
253  }
254  }
255 
256  for (int i = 0; i < idxs.size(); ++i) {
257  if (idxs[i] == xkey_) x = vals[i];
258  if (idxs[i] == ykey_) y = vals[i];
259  }
260  }
261 }
262 
263 PropertyPlot::PropertyPlot(const std::string& fluid_name, CoolProp::parameters ykey, CoolProp::parameters xkey, TPLimits tp_limits)
264  : xkey_(xkey),
265  ykey_(ykey) {
266  this->state_ = Detail::process_fluid_state(fluid_name);
267  this->critical_state_ = Detail::get_critical_point(state_);
268 
269  xaxis.scale = Detail::default_scale(xkey);
270  yaxis.scale = Detail::default_scale(ykey);
271 
272  // We are just assuming that all inputs and outputs are in SI units. We
273  // take care of any conversions before calling the library and after
274  // getting the results.
275  int out1, out2;
276  axis_pair_ = CoolProp::generate_update_pair(xkey, 0, ykey, 1, out1, out2);
277  swap_axis_inputs_for_update_ = (out1 == 1);
278 
279  const double HI_FACTOR = 2.25; // Upper default limits: HI_FACTOR*T_crit and HI_FACTOR*p_crit
280  const double LO_FACTOR = 1.01; // Lower default limits: LO_FACTOR*T_triple and LO_FACTOR*p_triple
281  switch (tp_limits) {
282  case TPLimits::None: this->Tp_limits_ = {{Detail::NaN, Detail::NaN}, {Detail::NaN, Detail::NaN}}; break;
283  case TPLimits::Def: this->Tp_limits_ = {{LO_FACTOR, HI_FACTOR}, {LO_FACTOR, HI_FACTOR}}; break;
284  case TPLimits::Achp: this->Tp_limits_ = {{173.15, 493.15}, {0.25e5, HI_FACTOR}}; break;
285  case TPLimits::Orc: this->Tp_limits_ = {{273.15, 673.15}, {0.25e5, HI_FACTOR}}; break;
286  }
287 
288  Range2D ranges = get_axis_limits();
289  xaxis.range = ranges.x;
290  yaxis.range = ranges.y;
291 }
292 
293 Range PropertyPlot::isoline_range(CoolProp::parameters key) const {
294  if (key == CoolProp::iQ)
295  return {0, 1};
296  else
297  return get_axis_limits(key, CoolProp::iT).x;
298 }
299 
300 Isolines PropertyPlot::calc_isolines(CoolProp::parameters key, const std::vector<double>& values, int points) const {
301  std::vector<double> xvals = generate_values_in_range(xaxis.scale, xaxis.range, points);
302  std::vector<double> yvals = generate_values_in_range(yaxis.scale, yaxis.range, points);
303 
304  Isolines lines;
305  for (double val : values) {
306  Isoline line(key, xkey_, ykey_, val, state_);
307  line.calc_range(xvals, yvals);
308  lines.push_back(line);
309  }
310  return lines;
311 }
312 
313 std::vector<CoolProp::parameters> PropertyPlot::supported_isoline_keys() const {
314  // taken from PropertyPlot::calc_isolines when called with iso_type='all'
315  std::vector<CoolProp::parameters> keys;
316  for (auto it = Detail::xy_switch.begin(); it != Detail::xy_switch.end(); ++it) {
317  const std::map<int, Detail::IsolineSupported>& supported = it->second;
318  auto supported_xy = supported.find(ykey_ * 10 + xkey_);
319  if (supported_xy != supported.end() && supported_xy->second != Detail::IsolineSupported::No)
320  keys.push_back(it->first);
321  }
322  return keys;
323 }
324 
325 double PropertyPlot::value_at(CoolProp::parameters key, double xvalue, double yvalue, CoolProp::phases phase) const {
326  if (key == xkey_) return xvalue;
327  if (key == ykey_) return yvalue;
328 
329  try {
330  if (swap_axis_inputs_for_update_)
331  std::swap(xvalue, yvalue);
332  state_->specify_phase(phase);
333  state_->update(axis_pair_, xvalue, yvalue);
334  switch (key) {
335  case CoolProp::iT: return state_->T();
336  case CoolProp::iP: return state_->p();
337  case CoolProp::iDmass: return state_->rhomass();
338  case CoolProp::iHmass: return state_->hmass();
339  case CoolProp::iSmass: return state_->smass();
340  case CoolProp::iUmass: return state_->umass();
341  case CoolProp::iQ: return state_->Q();
342  default: return Detail::NaN;
343  }
344  } catch (...) {
345  return Detail::NaN;
346  }
347 }
348 
349 Range PropertyPlot::get_sat_bounds(CoolProp::parameters key) const {
350  double s = 1e-7;
351  double t_small = critical_state_->keyed_output(CoolProp::iT) * s;
352  double p_small = critical_state_->keyed_output(CoolProp::iP) * s;
353 
354  double t_triple = state_->trivial_keyed_output(CoolProp::iT_triple);
355  double t_min;
356  try {
357  t_min = state_->trivial_keyed_output(CoolProp::iT_min);
358  } catch (...) {
359  t_min = t_triple;
360  }
361  state_->update(CoolProp::QT_INPUTS, 0, std::max(t_triple, t_min) + t_small);
362  if (key == CoolProp::iP)
363  return {state_->keyed_output(CoolProp::iP) + p_small, critical_state_->keyed_output(CoolProp::iP) - p_small};
364  else if (key == CoolProp::iT)
365  return {state_->keyed_output(CoolProp::iT) + t_small, critical_state_->keyed_output(CoolProp::iT) - t_small};
366  else
367  throw CoolProp::ValueError("Invalid key");
368 }
369 
370 PropertyPlot::Range2D PropertyPlot::get_Tp_limits() const {
371  Range t = Tp_limits_.T;
372  Range p = Tp_limits_.p;
373  Range tsat = get_sat_bounds(CoolProp::iT);
374  Range psat = get_sat_bounds(CoolProp::iP);
375 
376  const double ID_FACTOR = 10.0; // Values below this number are interpreted as factors
377  if (std::isnan(t.min)) t.min = 0.0;
378  else if (t.min < ID_FACTOR) t.min *= tsat.min;
379  if (std::isnan(t.max)) t.max = 1e6;
380  else if (t.max < ID_FACTOR) t.max *= tsat.max;
381  if (std::isnan(p.min)) p.min = 0.0;
382  else if (p.min < ID_FACTOR) p.min *= psat.min;
383  if (std::isnan(p.max)) p.max = 1e10;
384  else if (p.max < ID_FACTOR) p.max *= psat.max;
385 
386  try { t.min = std::max(t.min, state_->trivial_keyed_output(CoolProp::iT_min)); } catch (...) {}
387  try { t.max = std::min(t.max, state_->trivial_keyed_output(CoolProp::iT_max)); } catch (...) {}
388  try { p.min = std::max(p.min, state_->trivial_keyed_output(CoolProp::iP_min)); } catch (...) {}
389  try { p.max = std::min(p.max, state_->trivial_keyed_output(CoolProp::iP_max)); } catch (...) {}
390  return {t, p};
391 }
392 
393 PropertyPlot::Range2D PropertyPlot::get_axis_limits(CoolProp::parameters xkey, CoolProp::parameters ykey, bool autoscale) const {
394  if (xkey == CoolProp::parameters::iundefined_parameter) xkey = this->xkey_;
395  if (ykey == CoolProp::parameters::iundefined_parameter) ykey = this->ykey_;
396 
397  if (xkey != this->xkey_ || ykey != this->ykey_ || autoscale) {
398  Range2D tp_limits = get_Tp_limits();
399  Range xrange = {std::numeric_limits<double>::max(), std::numeric_limits<double>::lowest()};
400  Range yrange = {std::numeric_limits<double>::max(), std::numeric_limits<double>::lowest()};
401 
402  for (double T : {tp_limits.T.min, tp_limits.T.max}) {
403  for (double p : {tp_limits.p.min, tp_limits.p.max}) {
404  try {
405  state_->update(CoolProp::PT_INPUTS, p, T);
406  double x = state_->keyed_output(xkey);
407  double y = state_->keyed_output(ykey);
408  xrange.min = std::min(xrange.min, x);
409  xrange.max = std::max(xrange.max, x);
410  yrange.min = std::min(yrange.min, y);
411  yrange.max = std::max(yrange.max, y);
412  } catch (...) { }
413  }
414  }
415  return {xrange, yrange};
416  } else {
417  return {xaxis.range, yaxis.range};
418  }
419 }
420 
421 } /* namespace Plot */
422 } /* namespace CoolProp */
423 
424 
425 #ifdef ENABLE_CATCH
426 # include <catch2/catch_all.hpp>
427 # include <catch2/matchers/catch_matchers_floating_point.hpp>
428 
429 using Catch::Matchers::WithinAbs;
430 using Catch::Matchers::WithinRel;
431 
432 TEST_CASE("Check value_at for p-h plots", "[Plot]") {
433  CoolProp::Plot::PropertyPlot plot("R134a", CoolProp::iP, CoolProp::iHmass, CoolProp::Plot::TPLimits::Achp);
434 
435  CHECK_THAT(plot.value_at(CoolProp::iP, 300000/*Pa*/, 200000/*J/kg*/), WithinAbs(200000, 1e-10));
436  CHECK_THAT(plot.value_at(CoolProp::iHmass, 300000, 200000), WithinAbs(300000, 1e-10));
437  CHECK_THAT(plot.value_at(CoolProp::iT, 300000, 200000), WithinAbs(263.07372753976694, 1e-10));
438  CHECK_THAT(plot.value_at(CoolProp::iQ, 300000, 200000), WithinAbs(0.55044347874344737, 1e-10));
439 }
440 
441 TEST_CASE("Check that the isolines are the same as from Python", "[Plot]") {
442  CoolProp::Plot::PropertyPlot plot("HEOS::R134a", CoolProp::iP, CoolProp::iHmass, CoolProp::Plot::TPLimits::Achp);
443  const int isoline_count = 5;
444  const int points_per_isoline = 5;
445 
446  // CHECK(plot.xkey_ == CoolProp::iHmass);
447  // CHECK(plot.ykey_ == CoolProp::iP);
448 
449  CHECK(plot.xaxis.scale == CoolProp::Plot::Scale::Lin);
450  CHECK(plot.yaxis.scale == CoolProp::Plot::Scale::Log);
451  CHECK_THAT(plot.xaxis.min, WithinAbs(75373.1, 1));
452  CHECK_THAT(plot.xaxis.max, WithinAbs(577605, 1));
453  CHECK_THAT(plot.yaxis.min, WithinAbs(25000, 1));
454  CHECK_THAT(plot.yaxis.max, WithinAbs(9.133e6, 1));
455 
456  std::vector<CoolProp::parameters> iso_types = plot.supported_isoline_keys();
457  REQUIRE(iso_types.size() == 4);
458  CHECK(iso_types[0] == CoolProp::iT);
459  CHECK(iso_types[1] == CoolProp::iQ);
460  CHECK(iso_types[2] == CoolProp::iDmass);
461  CHECK(iso_types[3] == CoolProp::iSmass);
462 
463  {
464  // Q isolines
465  CoolProp::Plot::Range q_range = plot.isoline_range(CoolProp::iQ);
466  std::vector<double> q_values = CoolProp::Plot::generate_values_in_range(CoolProp::iQ, q_range, isoline_count);
467  CoolProp::Plot::Isolines q_isolines = plot.calc_isolines(CoolProp::iQ, q_values, points_per_isoline);
468  REQUIRE(q_isolines.size() == isoline_count);
469  CHECK_THAT(q_isolines[0].value, WithinAbs(0.0, 1e-10));
470  CHECK_THAT(q_isolines[1].value, WithinAbs(0.25, 1e-10));
471  CHECK_THAT(q_isolines[2].value, WithinAbs(0.5, 1e-10));
472  CHECK_THAT(q_isolines[3].value, WithinAbs(0.75, 1e-10));
473  CHECK_THAT(q_isolines[4].value, WithinAbs(1.0, 1e-10));
474  const double expected_x[isoline_count][points_per_isoline] = {
475  {71455.08256999, 132939.99472497, 198497.0525314, 271576.58908888, 389440.73899019},
476  {137326.83116781, 191267.05172013, 248359.90642508, 309538.95484829, 389511.40516519},
477  {203198.57976564, 249594.1087153, 298222.76031877, 347501.3206077, 389582.0713402},
478  {269070.32836347, 307921.16571046, 348085.61421246, 385463.68636711, 389652.73751521},
479  {334942.07696129, 366248.22270562, 397948.46810615, 423426.05212652, 389723.40369022}
480  };
481  const double expected_y[isoline_count][points_per_isoline] = {
482  {3.89567060e+02, 2.58505756e+04, 2.81105747e+05, 1.31691170e+06, 4.05910826e+06},
483  {3.89567060e+02, 2.58505756e+04, 2.81105747e+05, 1.31691170e+06, 4.05910826e+06},
484  {3.89567060e+02, 2.58505756e+04, 2.81105747e+05, 1.31691170e+06, 4.05910826e+06},
485  {3.89567060e+02, 2.58505756e+04, 2.81105747e+05, 1.31691170e+06, 4.05910826e+06},
486  {3.89567060e+02, 2.58505756e+04, 2.81105747e+05, 1.31691170e+06, 4.05910826e+06}
487  };
488  for (int i = 0; i < q_isolines.size(); ++i) {
489  REQUIRE(q_isolines[i].size() == points_per_isoline);
490  for (int j = 0; j < q_isolines[i].size(); ++j) {
491  CHECK_THAT(q_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
492  CHECK_THAT(q_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
493  }
494  }
495  }
496  {
497  // T isolines
498  CoolProp::Plot::Range t_range = plot.isoline_range(CoolProp::iT);
499  std::vector<double> t_values = CoolProp::Plot::generate_values_in_range(CoolProp::iT, t_range, isoline_count);
500  CoolProp::Plot::Isolines t_isolines = plot.calc_isolines(CoolProp::iT, t_values, points_per_isoline);
501  REQUIRE(t_isolines.size() == isoline_count);
502  CHECK_THAT(t_isolines[0].value, WithinAbs(173.15, 1e-10));
503  CHECK_THAT(t_isolines[1].value, WithinAbs(243.6125, 1e-10));
504  CHECK_THAT(t_isolines[2].value, WithinAbs(314.075, 1e-10));
505  CHECK_THAT(t_isolines[3].value, WithinAbs(384.5375, 1e-10));
506  CHECK_THAT(t_isolines[4].value, WithinAbs(455.0, 1e-10));
507  const double expected_x[isoline_count][points_per_isoline] = {
508  {75373.12689908, 75410.99061368, 75576.57734102, 76301.46320034, 79487.71936123},
509  {382785.23058756, 161389.44197265, 161516.21527543, 162076.96181603, 164636.92352411},
510  {439466.64984328, 438148.19040202, 431912.24266074, 257605.32897193, 257512.80690587},
511  {504550.62606561, 503783.53950874, 500331.68636519, 482707.98489058, 366959.25257106},
512  {577604.59497521, 577097.07174302, 574850.21206939, 564444.22079795, 507878.75380996}
513  };
514  const double expected_y[isoline_count][points_per_isoline] = {
515  {25000.0, 109297.03270763, 477833.65434772, 2089032.02192197, 9133000.04909175},
516  {25000.0, 109297.03270763, 477833.65434772, 2089032.02192197, 9133000.04909175},
517  {25000.0, 109297.03270763, 477833.65434772, 2089032.02192197, 9133000.04909175},
518  {25000.0, 109297.03270763, 477833.65434772, 2089032.02192197, 9133000.04909175},
519  {25000.0, 109297.03270763, 477833.65434772, 2089032.02192197, 9133000.04909175}
520  };
521  for (int i = 0; i < t_isolines.size(); ++i) {
522  REQUIRE(t_isolines[i].size() == points_per_isoline);
523  for (int j = 0; j < t_isolines[i].size(); ++j) {
524  CHECK_THAT(t_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
525  CHECK_THAT(t_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
526  }
527  }
528  }
529  {
530  // S isolines
531  CoolProp::Plot::Range s_range = plot.isoline_range(CoolProp::iSmass);
532  std::vector<double> s_values = CoolProp::Plot::generate_values_in_range(CoolProp::iSmass, s_range, isoline_count);
533  CoolProp::Plot::Isolines s_isolines = plot.calc_isolines(CoolProp::iSmass, s_values, points_per_isoline);
534  REQUIRE(s_isolines.size() == isoline_count);
535  CHECK_THAT(s_isolines[0].value, WithinAbs(426.0098553415565, 1e-10));
536  CHECK_THAT(s_isolines[1].value, WithinAbs(925.2753143522199, 1e-10));
537  CHECK_THAT(s_isolines[2].value, WithinAbs(1424.540773362883, 1e-10));
538  CHECK_THAT(s_isolines[3].value, WithinAbs(1923.8062323735467, 1e-10));
539  CHECK_THAT(s_isolines[4].value, WithinAbs(2423.07169138421, 1e-10));
540  const double expected_x[isoline_count][points_per_isoline] = {
541  {73758.20064883, 73811.34928194, 74043.68191583, 75058.90103546, 79487.71936135},
542  {176257.41124603, 179794.86552304, 180290.38376303, 181487.99233203, 186690.75978367},
543  {286286.21675221, 303984.6485323, 321692.18498473, 335551.56116185, 344087.52152244},
544  {399372.58517377, 433400.13264058, 471964.37092222, 513835.0906295, 555823.55477128},
545  {577604.59497522, 635257.81956317, 698998.52697158, 768744.13132575, std::nan("")}
546  };
547  const double expected_y[isoline_count][points_per_isoline] = {
548  {25000.0, 109297.03270763, 477833.65434772, 2089032.02192197, 9133000.04909175},
549  {25000.0, 109297.03270763, 477833.65434772, 2089032.02192197, 9133000.04909175},
550  {25000.0, 109297.03270763, 477833.65434772, 2089032.02192197, 9133000.04909175},
551  {25000.0, 109297.03270763, 477833.65434772, 2089032.02192197, 9133000.04909175},
552  {25000.0, 109297.03270763, 477833.65434772, 2089032.02192197, 9133000.04909175}
553  };
554  for (int i = 0; i < s_isolines.size(); ++i) {
555  REQUIRE(s_isolines[i].size() == points_per_isoline);
556  for (int j = 0; j < s_isolines[i].size(); ++j) {
557  if (std::isnan(s_isolines[i].x[j]))
558  CHECK(std::isnan(expected_x[i][j]));
559  else
560  CHECK_THAT(s_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
561  CHECK_THAT(s_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
562  }
563  }
564  }
565  {
566  // D isolines
567  CoolProp::Plot::Range d_range = plot.isoline_range(CoolProp::iDmass);
568  std::vector<double> d_values = CoolProp::Plot::generate_values_in_range(CoolProp::iDmass, d_range, isoline_count);
569  CoolProp::Plot::Isolines d_isolines = plot.calc_isolines(CoolProp::iDmass, d_values, points_per_isoline);
570  REQUIRE(d_isolines.size() == isoline_count);
571  CHECK_THAT(d_isolines[0].value, WithinAbs(0.6749779869915704, 1e-10));
572  CHECK_THAT(d_isolines[1].value, WithinAbs(4.704765330619733, 1e-10));
573  CHECK_THAT(d_isolines[2].value, WithinAbs(32.793390662794806, 1e-10));
574  CHECK_THAT(d_isolines[3].value, WithinAbs(228.57813208316188, 1e-10));
575  CHECK_THAT(d_isolines[4].value, WithinAbs(1593.2467308391022, 1e-10));
576  const double expected_x[isoline_count][points_per_isoline] = {
577  {5.77604595e+05, 3.65397965e+06, 3.84283606e+07, 4.40954815e+08, 5.22494051e+09},
578  {2.02365849e+05, 4.19227802e+05, 1.84512838e+06, 1.78590373e+07, 2.01674048e+08},
579  {1.42114492e+05, 2.04387395e+05, 3.51213643e+05, 1.00846111e+06, 8.12669299e+06},
580  {1.33470419e+05, 1.72415441e+05, 2.35381904e+05, 3.57488786e+05, 6.69475618e+05},
581  {7.05185202e+04, 7.06013993e+04, 7.09637630e+04, 7.25484894e+04, 7.94877194e+04}
582  };
583  const double expected_y[isoline_count][points_per_isoline] = {
584  {25000.0, 109297.03270763, 477833.65434772, 2089032.02192197, 9133000.04909175},
585  {25000.0, 109297.03270763, 477833.65434772, 2089032.02192197, 9133000.04909175},
586  {25000.0, 109297.03270763, 477833.65434772, 2089032.02192197, 9133000.04909175},
587  {25000.0, 109297.03270763, 477833.65434772, 2089032.02192197, 9133000.04909175},
588  {25000.0, 109297.03270763, 477833.65434772, 2089032.02192197, 9133000.04909175}
589  };
590  for (int i = 0; i < d_isolines.size(); ++i) {
591  REQUIRE(d_isolines[i].size() == points_per_isoline);
592  for (int j = 0; j < d_isolines[i].size(); ++j) {
593  CHECK_THAT(d_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
594  CHECK_THAT(d_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
595  }
596  }
597  }
598 }
599 
600 TEST_CASE("Basic TS Plot has same output as Python", "[Plot]") {
601  CoolProp::Plot::PropertyPlot plot("HEOS::R134a", CoolProp::iT, CoolProp::iSmass, CoolProp::Plot::TPLimits::Achp);
602  const int isoline_count = 5;
603  const int points_per_isoline = 5;
604 
605  // CHECK(plot.xkey_ == CoolProp::iSmass);
606  // CHECK(plot.ykey_ == CoolProp::iT);
607 
608  CHECK(plot.xaxis.scale == CoolProp::Plot::Scale::Lin);
609  CHECK(plot.yaxis.scale == CoolProp::Plot::Scale::Lin);
610  CHECK_THAT(plot.xaxis.min, WithinAbs(426, 1));
611  CHECK_THAT(plot.xaxis.max, WithinAbs(2423, 1));
612  CHECK_THAT(plot.yaxis.min, WithinAbs(173, 1));
613  CHECK_THAT(plot.yaxis.max, WithinAbs(455, 1));
614 
615  std::vector<CoolProp::parameters> iso_types = plot.supported_isoline_keys();
616  REQUIRE(iso_types.size() == 4);
617  CHECK(iso_types[0] == CoolProp::iP);
618  CHECK(iso_types[1] == CoolProp::iQ);
619  CHECK(iso_types[2] == CoolProp::iDmass);
620  CHECK(iso_types[3] == CoolProp::iHmass);
621 
622  {
623  // Q isolines
624  CoolProp::Plot::Range q_range = plot.isoline_range(CoolProp::iQ);
625  std::vector<double> q_values = CoolProp::Plot::generate_values_in_range(CoolProp::iQ, q_range, isoline_count);
626  CoolProp::Plot::Isolines q_isolines = plot.calc_isolines(CoolProp::iQ, q_values, points_per_isoline);
627  REQUIRE(q_isolines.size() == isoline_count);
628  CHECK_THAT(q_isolines[0].value, WithinAbs(0.0, 1e-10));
629  CHECK_THAT(q_isolines[1].value, WithinAbs(0.25, 1e-10));
630  CHECK_THAT(q_isolines[2].value, WithinAbs(0.5, 1e-10));
631  CHECK_THAT(q_isolines[3].value, WithinAbs(0.75, 1e-10));
632  CHECK_THAT(q_isolines[4].value, WithinAbs(1.0, 1e-10));
633 
634  const double expected_x[isoline_count][points_per_isoline] = {
635  {412.61753823, 728.71208313, 994.51958846, 1237.3122956, 1561.56967158},
636  {800.44043827, 992.70703491, 1177.81867482, 1354.79919476, 1561.75851162},
637  {1188.26333832, 1256.70198669, 1361.11776119, 1472.28609393, 1561.94735166},
638  {1576.08623836, 1520.69693847, 1544.41684755, 1589.77299309, 1562.1361917},
639  {1963.9091384, 1784.69189025, 1727.71593392, 1707.25989226, 1562.32503174}
640  };
641  const double expected_y[isoline_count][points_per_isoline] = {
642  {169.85007484, 220.94004678, 272.03001871, 323.11999064, 374.20996258},
643  {169.85007484, 220.94004678, 272.03001871, 323.11999064, 374.20996258},
644  {169.85007484, 220.94004678, 272.03001871, 323.11999064, 374.20996258},
645  {169.85007484, 220.94004678, 272.03001871, 323.11999064, 374.20996258},
646  {169.85007484, 220.94004678, 272.03001871, 323.11999064, 374.20996258}
647  };
648 
649  for(int i = 0; i < q_isolines.size(); i++) {
650  REQUIRE(q_isolines[i].size() == points_per_isoline);
651  for(int j = 0; j < q_isolines[i].size(); j++) {
652  CHECK_THAT(q_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
653  CHECK_THAT(q_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
654  }
655  }
656  }
657  {
658  // P isolines
659  CoolProp::Plot::Range p_range = plot.isoline_range(CoolProp::iP);
660  std::vector<double> p_values = CoolProp::Plot::generate_values_in_range(CoolProp::iP, p_range, isoline_count);
661  CoolProp::Plot::Isolines p_isolines = plot.calc_isolines(CoolProp::iP, p_values, points_per_isoline);
662  REQUIRE(p_isolines.size() == isoline_count);
663  CHECK_THAT(p_isolines[0].value, WithinAbs(25000.000000000007, 1e-8));
664  CHECK_THAT(p_isolines[1].value, WithinAbs(109297.03270763098, 1e-8));
665  CHECK_THAT(p_isolines[2].value, WithinAbs(477833.65434771683, 1e-8));
666  CHECK_THAT(p_isolines[3].value, WithinAbs(2089032.0219219688, 1e-8));
667  CHECK_THAT(p_isolines[4].value, WithinAbs(9133000.049091753, 1e-8));
668 
669  const double expected_x[isoline_count][points_per_isoline] = {
670  {426.00985534, 925.27531435, 1424.54077336, 1923.80623237, 2423.07169138},
671  {426.00985534, 925.27531435, 1424.54077336, 1923.80623237, 2423.07169138},
672  {426.00985534, 925.27531435, 1424.54077336, 1923.80623237, 2423.07169138},
673  {426.00985534, 925.27531435, 1424.54077336, 1923.80623237, 2423.07169138},
674  {426.00985534, 925.27531435, 1424.54077336, 1923.80623237, 2423.07169138}
675  };
676  const double expected_y[isoline_count][points_per_isoline] = {
677  {171.78612656, 220.38136931, 220.38136931, 265.36250878, 455.0},
678  {171.7989644, 248.7449928, 248.7449928, 308.63364605, 506.38739121},
679  {171.85504128, 258.19575097, 287.47240852, 355.96420961, 560.29233808},
680  {172.09927987, 258.74250558, 342.56046108, 411.32270988, 618.0350615},
681  {173.15, 261.02100275, 371.32673696, 484.42563591, std::nan("")}
682  };
683 
684  for(int i = 0; i < p_isolines.size(); i++) {
685  REQUIRE(p_isolines[i].size() == points_per_isoline);
686  for(int j = 0; j < p_isolines[i].size(); j++) {
687  if (std::isnan(expected_y[i][j])) {
688  CHECK(std::isnan(p_isolines[i].y[j]));
689  } else {
690  CHECK_THAT(p_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
691  CHECK_THAT(p_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
692  }
693  }
694  }
695  }
696  {
697  // H isolines
698  CoolProp::Plot::Range h_range = plot.isoline_range(CoolProp::iHmass);
699  std::vector<double> h_values = CoolProp::Plot::generate_values_in_range(CoolProp::iHmass, h_range, isoline_count);
700  CoolProp::Plot::Isolines h_isolines = plot.calc_isolines(CoolProp::iHmass, h_values, points_per_isoline);
701  REQUIRE(h_isolines.size() == isoline_count);
702  CHECK_THAT(h_isolines[0].value, WithinAbs(75373.12689908482, 1e-10));
703  CHECK_THAT(h_isolines[1].value, WithinAbs(200930.99391811725, 1e-10));
704  CHECK_THAT(h_isolines[2].value, WithinAbs(326488.8609371497, 1e-10));
705  CHECK_THAT(h_isolines[3].value, WithinAbs(452046.72795618215, 1e-10));
706  CHECK_THAT(h_isolines[4].value, WithinAbs(577604.5949752146, 1e-10));
707 
708  const double expected_x[isoline_count][points_per_isoline] = {
709  {426.00985534, 925.27531435, 1424.54077336, 1923.80623237, 2423.07169138},
710  {426.00985534, 925.27531435, 1424.54077336, 1923.80623237, 2423.07169138},
711  {426.00985534, 925.27531435, 1424.54077336, 1923.80623237, 2423.07169138},
712  {426.00985534, 925.27531435, 1424.54077336, 1923.80623237, 2423.07169138},
713  {426.00985534, 925.27531435, 1424.54077336, 1923.80623237, 2423.07169138}
714  };
715  const double expected_y[isoline_count][points_per_isoline] = {
716  {172.17461409, std::nan(""), std::nan(""), std::nan(""), std::nan("")},
717  {196.07460266, 266.63119128, std::nan(""), std::nan(""), std::nan("")},
718  {213.66474036, 299.9847036, 301.72638717, std::nan(""), std::nan("")},
719  {228.4112647, 322.84362137, 426.78791645, 331.52116588, 328.04216753},
720  {241.56832462, 341.66140042, 458.59389239, std::nan(""), 455.0}
721  };
722 
723  for(int i = 0; i < h_isolines.size(); i++) {
724  REQUIRE(h_isolines[i].size() == points_per_isoline);
725  for(int j = 0; j < h_isolines[i].size(); j++) {
726  if (std::isnan(expected_y[i][j])) {
727  CHECK(std::isnan(h_isolines[i].y[j]));
728  } else {
729  CHECK_THAT(h_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
730  CHECK_THAT(h_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
731  }
732  }
733  }
734  }
735  {
736  // D isolines
737  CoolProp::Plot::Range d_range = plot.isoline_range(CoolProp::iDmass);
738  std::vector<double> d_values = CoolProp::Plot::generate_values_in_range(CoolProp::iDmass, d_range, isoline_count);
739  CoolProp::Plot::Isolines d_isolines = plot.calc_isolines(CoolProp::iDmass, d_values, points_per_isoline);
740  REQUIRE(d_isolines.size() == isoline_count);
741  CHECK_THAT(d_isolines[0].value, WithinAbs(0.6749779869915704, 1e-10));
742  CHECK_THAT(d_isolines[1].value, WithinAbs(4.704765330619733, 1e-10));
743  CHECK_THAT(d_isolines[2].value, WithinAbs(32.793390662794806, 1e-10));
744  CHECK_THAT(d_isolines[3].value, WithinAbs(228.57813208316188, 1e-10));
745  CHECK_THAT(d_isolines[4].value, WithinAbs(1593.2467308391022, 1e-10));
746 
747  const double expected_x[isoline_count][points_per_isoline] = {
748  {524.17387831, 1911.09303198, 2092.95299736, 2262.71394473, 2423.07169138},
749  {448.10309047, 1715.11962047, 1932.46628376, 2103.15612883, 2263.90954344},
750  {437.18945214, 972.48976631, 1758.36242394, 1935.75229847, 2099.20644384},
751  {435.62370654, 865.94698069, 1292.02342121, 1720.27748872, 1899.3816054},
752  {426.00985534, 710.87738683, 946.96900373, 1151.9181105, 1335.56535146}
753  };
754  const double expected_y[isoline_count][points_per_isoline] = {
755  {173.15, 243.6125, 314.075, 384.5375, 455.0},
756  {173.15, 243.6125, 314.075, 384.5375, 455.0},
757  {173.15, 243.6125, 314.075, 384.5375, 455.0},
758  {173.15, 243.6125, 314.075, 384.5375, 455.0},
759  {173.15, 243.6125, 314.075, 384.5375, 455.0}
760  };
761 
762  for(int i = 0; i < d_isolines.size(); i++) {
763  REQUIRE(d_isolines[i].size() == points_per_isoline);
764  for(int j = 0; j < d_isolines[i].size(); j++) {
765  CHECK_THAT(d_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
766  CHECK_THAT(d_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
767  }
768  }
769  }
770 }
771 
772 #endif /* ENABLE_CATCH */
Generated on Mon Mar 17 2025 04:25:23 for CoolProp by doxygen 1.9.1