CoolProp  6.6.1dev
An open-source fluid property and humid air property database
TabularBackends.cpp
Go to the documentation of this file.
1 
2 #if !defined(NO_TABULAR_BACKENDS)
3 
4 # include "TabularBackends.h"
5 # include "CoolProp.h"
6 # include <sstream>
7 # include "time.h"
8 # include "miniz.h"
9 # include <fstream>
10 
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);
22 
23 static CoolProp::TabularDataLibrary library;
24 
25 namespace CoolProp {
26 
33 template <typename T>
34 void load_table(T& table, const std::string& path_to_tables, const std::string& filename) {
35 
36  double tic = clock();
37  std::string path_to_table = path_to_tables + "/" + filename;
38  if (get_debug_level() > 0) {
39  std::cout << format("Loading table: %s", path_to_table.c_str()) << std::endl;
40  }
41  std::vector<char> raw;
42  try {
43  raw = get_binary_file_contents(path_to_table.c_str());
44  } catch (...) {
45  std::string err = format("Unable to load file %s", path_to_table.c_str());
46  if (get_debug_level() > 0) {
47  std::cout << "err:" << err << std::endl;
48  }
49  throw UnableToLoadError(err);
50  }
51  std::vector<unsigned char> newBuffer(raw.size() * 5);
52  uLong newBufferSize = static_cast<uLong>(newBuffer.size());
53  mz_ulong rawBufferSize = static_cast<mz_ulong>(raw.size());
54  int code;
55  do {
56  code = uncompress((unsigned char*)(&(newBuffer[0])), &newBufferSize, (unsigned char*)(&(raw[0])), rawBufferSize);
57  if (code == Z_BUF_ERROR) {
58  // Output buffer is too small, make it bigger and try again
59  newBuffer.resize(newBuffer.size() * 5);
60  newBufferSize = static_cast<uLong>(newBuffer.size());
61  } else if (code != 0) { // Something else, a big problem
62  std::string err = format("Unable to uncompress file %s with miniz code %d", path_to_table.c_str(), code);
63  if (get_debug_level() > 0) {
64  std::cout << "uncompress err:" << err << std::endl;
65  }
66  throw UnableToLoadError(err);
67  }
68  } while (code != 0);
69  // Copy the buffer from unsigned char to char (yuck)
70  std::vector<char> charbuffer(newBuffer.begin(), newBuffer.begin() + newBufferSize);
71  try {
72  msgpack::unpacked msg;
73  msgpack::unpack(msg, &(charbuffer[0]), charbuffer.size());
74  msgpack::object deserialized = msg.get();
75 
76  // Call the class' deserialize function; if it is an invalid table, it will cause an exception to be thrown
77  table.deserialize(deserialized);
78  double toc = clock();
79  if (get_debug_level() > 0) {
80  std::cout << format("Loaded table: %s in %g sec.", path_to_table.c_str(), (toc - tic) / CLOCKS_PER_SEC) << std::endl;
81  }
82  } catch (std::exception& e) {
83  std::string err = format("Unable to msgpack deserialize %s; err: %s", path_to_table.c_str(), e.what());
84  if (get_debug_level() > 0) {
85  std::cout << "err: " << err << std::endl;
86  }
87  throw UnableToLoadError(err);
88  }
89 }
90 template <typename T>
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);
101  ofs2.close();
102 
103  if (CoolProp::get_config_bool(SAVE_RAW_TABLES)) {
104  std::ofstream ofs(tabPath.c_str(), std::ofstream::binary);
105  ofs.write(sbuf.data(), sbuf.size());
106  }
107 }
108 
109 } // namespace CoolProp
110 
111 void CoolProp::PureFluidSaturationTableData::build(shared_ptr<CoolProp::AbstractState>& AS) {
112  const bool debug = get_debug_level() > 5 || false;
113  if (debug) {
114  std::cout << format("***********************************************\n");
115  std::cout << format(" Saturation Table (%s) \n", AS->name().c_str());
116  std::cout << format("***********************************************\n");
117  }
118  resize(N);
119  // ------------------------
120  // Actually build the table
121  // ------------------------
122  CoolPropDbl Tmin = std::max(AS->Ttriple(), AS->Tmin());
123  AS->update(QT_INPUTS, 0, Tmin);
124  CoolPropDbl p_triple = AS->p();
125  CoolPropDbl p, pmin = p_triple, pmax = 0.9999 * AS->p_critical();
126  for (std::size_t i = 0; i < N - 1; ++i) {
127  if (i == 0) {
128  CoolProp::set_config_bool(DONT_CHECK_PROPERTY_LIMITS, true);
129  }
130  // Log spaced
131  p = exp(log(pmin) + (log(pmax) - log(pmin)) / (N - 1) * i);
132 
133  // Saturated liquid
134  try {
135  AS->update(PQ_INPUTS, p, 0);
136  pL[i] = p;
137  TL[i] = AS->T();
138  rhomolarL[i] = AS->rhomolar();
139  hmolarL[i] = AS->hmolar();
140  smolarL[i] = AS->smolar();
141  umolarL[i] = AS->umolar();
142  logpL[i] = log(p);
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) {
148  // That failed for some reason, go to the next pair
149  if (debug) {
150  std::cout << " " << e.what() << std::endl;
151  }
152  continue;
153  }
154  // Transport properties - if no transport properties, just keep going
155  try {
156  viscL[i] = AS->viscosity();
157  condL[i] = AS->conductivity();
158  logviscL[i] = log(viscL[i]);
159  } catch (std::exception& e) {
160  if (debug) {
161  std::cout << " " << e.what() << std::endl;
162  }
163  }
164  // Saturated vapor
165  try {
166  AS->update(PQ_INPUTS, p, 1);
167  pV[i] = p;
168  TV[i] = AS->T();
169  rhomolarV[i] = AS->rhomolar();
170  hmolarV[i] = AS->hmolar();
171  smolarV[i] = AS->smolar();
172  umolarV[i] = AS->umolar();
173  logpV[i] = log(p);
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) {
179  // That failed for some reason, go to the next pair
180  if (debug) {
181  std::cout << " " << e.what() << std::endl;
182  }
183  continue;
184  }
185  // Transport properties - if no transport properties, just keep going
186  try {
187  viscV[i] = AS->viscosity();
188  condV[i] = AS->conductivity();
189  logviscV[i] = log(viscV[i]);
190  } catch (std::exception& e) {
191  if (debug) {
192  std::cout << " " << e.what() << std::endl;
193  }
194  }
195  if (i == 0) {
196  CoolProp::set_config_bool(DONT_CHECK_PROPERTY_LIMITS, false);
197  }
198  }
199  // Last point is at the critical point
200  AS->update(PQ_INPUTS, AS->p_critical(), 1);
201  std::size_t i = N - 1;
202 
203  pV[i] = AS->p();
204  TV[i] = AS->T();
205  rhomolarV[i] = AS->rhomolar();
206  hmolarV[i] = AS->hmolar();
207  smolarV[i] = AS->smolar();
208  umolarV[i] = AS->umolar();
209 
210  pL[i] = AS->p();
211  TL[i] = AS->T();
212  rhomolarL[i] = AS->rhomolar();
213  hmolarL[i] = AS->hmolar();
214  smolarL[i] = AS->smolar();
215  umolarL[i] = AS->umolar();
216 
217  logpV[i] = log(AS->p());
218  logrhomolarV[i] = log(rhomolarV[i]);
219 
220  logpL[i] = log(AS->p());
221  logrhomolarL[i] = log(rhomolarL[i]);
222 }
223 
224 void CoolProp::SinglePhaseGriddedTableData::build(shared_ptr<CoolProp::AbstractState>& AS) {
225  CoolPropDbl x, y;
226  const bool debug = get_debug_level() > 5 || false;
227 
228  resize(Nx, Ny);
229 
230  if (debug) {
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");
234  }
235  // ------------------------
236  // Actually build the table
237  // ------------------------
238  for (std::size_t i = 0; i < Nx; ++i) {
239  // Calculate the x value
240  if (logx) {
241  // Log spaced
242  x = exp(log(xmin) + (log(xmax) - log(xmin)) / (Nx - 1) * i);
243  } else {
244  // Linearly spaced
245  x = xmin + (xmax - xmin) / (Nx - 1) * i;
246  }
247  xvec[i] = x;
248  for (std::size_t j = 0; j < Ny; ++j) {
249  // Calculate the x value
250  if (logy) {
251  // Log spaced
252  y = exp(log(ymin) + (log(ymax / ymin)) / (Ny - 1) * j);
253  } else {
254  // Linearly spaced
255  y = ymin + (ymax - ymin) / (Ny - 1) * j;
256  }
257  yvec[j] = y;
258 
259  if (debug) {
260  std::cout << "x: " << x << " y: " << y << std::endl;
261  }
262 
263  // Generate the input pair
264  CoolPropDbl v1, v2;
265  input_pairs input_pair = generate_update_pair(xkey, x, ykey, y, v1, v2);
266 
267  // --------------------
268  // Update the state
269  // --------------------
270  try {
271  AS->update(input_pair, v1, v2);
272  if (!ValidNumber(AS->rhomolar())) {
273  throw ValueError("rhomolar is invalid");
274  }
275  } catch (std::exception& e) {
276  // That failed for some reason, go to the next pair
277  if (debug) {
278  std::cout << " " << e.what() << std::endl;
279  }
280  continue;
281  }
282 
283  // Skip two-phase states - they will remain as _HUGE holes in the table
284  if (is_in_closed_range(0.0, 1.0, AS->Q())) {
285  if (debug) {
286  std::cout << " 2Phase" << std::endl;
287  }
288  continue;
289  };
290 
291  // --------------------
292  // State variables
293  // --------------------
294  T[i][j] = AS->T();
295  p[i][j] = AS->p();
296  rhomolar[i][j] = AS->rhomolar();
297  hmolar[i][j] = AS->hmolar();
298  smolar[i][j] = AS->smolar();
299  umolar[i][j] = AS->umolar();
300 
301  // -------------------------
302  // Transport properties
303  // -------------------------
304  try {
305  visc[i][j] = AS->viscosity();
306  cond[i][j] = AS->conductivity();
307  } catch (std::exception&) {
308  // Failures will remain as holes in table
309  }
310 
311  // ----------------------------------------
312  // First derivatives of state variables
313  // ----------------------------------------
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);
326 
327  // ----------------------------------------
328  // Second derivatives of state variables
329  // ----------------------------------------
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);
348  }
349  }
350 }
351 std::string CoolProp::TabularBackend::path_to_tables(void) {
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]));
357  }
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;
362  }
363  return table_directory + AS->backend_name() + "(" + strjoin(components, "&") + ")";
364 }
365 
366 void CoolProp::TabularBackend::write_tables() {
367  std::string path_to_tables = this->path_to_tables();
368  make_dirs(path_to_tables);
369  bool loaded = false;
370  dataset = library.get_set_of_tables(this->AS, loaded);
371  PackablePhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
372  PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
373  SinglePhaseGriddedTableData& single_phase_logph = dataset->single_phase_logph;
374  SinglePhaseGriddedTableData& single_phase_logpT = dataset->single_phase_logpT;
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");
379 }
380 void CoolProp::TabularBackend::load_tables() {
381  bool loaded = false;
382  dataset = library.get_set_of_tables(this->AS, loaded);
383  if (loaded == false) {
384  throw UnableToLoadError("Could not load tables");
385  }
386  if (get_debug_level() > 0) {
387  std::cout << "Tables loaded" << std::endl;
388  }
389 }
390 
391 CoolPropDbl CoolProp::TabularBackend::calc_saturated_vapor_keyed_output(parameters key) {
392  PhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
393  PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
394  double factor = 1.0;
395  mass_to_molar(key, factor, molar_mass());
396  if (is_mixture) {
397  return phase_envelope_sat(phase_envelope, key, iP, _p) * factor;
398  } else {
399  return pure_saturation.evaluate(key, _p, 1, cached_saturation_iL, cached_saturation_iV) * factor;
400  }
401 }
402 CoolPropDbl CoolProp::TabularBackend::calc_saturated_liquid_keyed_output(parameters key) {
403  PhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
404  PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
405  double factor = 1.0;
406  mass_to_molar(key, factor, molar_mass());
407  if (is_mixture) {
408  return phase_envelope_sat(phase_envelope, key, iP, _p) * factor;
409  } else {
410  return pure_saturation.evaluate(key, _p, 0, cached_saturation_iL, cached_saturation_iV) * factor;
411  }
412 };
413 
414 CoolPropDbl CoolProp::TabularBackend::calc_p(void) {
415  PhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
416  if (using_single_phase_table) {
417  return _p;
418  } else {
419  if (is_mixture) {
420  return phase_envelope_sat(phase_envelope, iP, iT, _T);
421  } else {
422  return _p;
423  }
424  }
425 }
426 CoolPropDbl CoolProp::TabularBackend::calc_T(void) {
427  PhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
428  PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
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:
434  return _T;
435  case SELECTED_NO_TABLE:
436  throw ValueError("table not selected");
437  }
438  return _HUGE; // not needed, will never be hit, just to make compiler happy
439  } else {
440  if (is_mixture) {
441  return phase_envelope_sat(phase_envelope, iT, iP, _p);
442  } else {
443  if (ValidNumber(_T)) {
444  return _T;
445  } else {
446  return pure_saturation.evaluate(iT, _p, _Q, cached_saturation_iL, cached_saturation_iV);
447  }
448  }
449  }
450 }
451 CoolPropDbl CoolProp::TabularBackend::calc_rhomolar(void) {
452  PhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
453  PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
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:
461  throw ValueError("table not selected");
462  }
463  return _HUGE; // not needed, will never be hit, just to make compiler happy
464  } else {
465  if (is_mixture) {
466  return phase_envelope_sat(phase_envelope, iDmolar, iP, _p);
467  } else {
468  return pure_saturation.evaluate(iDmolar, _p, _Q, cached_saturation_iL, cached_saturation_iV);
469  }
470  }
471 }
472 CoolPropDbl CoolProp::TabularBackend::calc_hmolar(void) {
473  PhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
474  PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
475  if (using_single_phase_table) {
476  switch (selected_table) {
477  case SELECTED_PH_TABLE:
478  return _hmolar;
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:
482  throw ValueError("table not selected");
483  }
484  return _HUGE; // not needed, will never be hit, just to make compiler happy
485  } else {
486  if (is_mixture) {
487  return phase_envelope_sat(phase_envelope, iHmolar, iP, _p);
488  } else {
489  return pure_saturation.evaluate(iHmolar, _p, _Q, cached_saturation_iL, cached_saturation_iV);
490  }
491  }
492 }
493 CoolPropDbl CoolProp::TabularBackend::calc_smolar(void) {
494  PhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
495  PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
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:
503  throw ValueError("table not selected");
504  }
505  return _HUGE; // not needed, will never be hit, just to make compiler happy
506  } else {
507  if (is_mixture) {
508  return phase_envelope_sat(phase_envelope, iSmolar, iP, _p);
509  } else {
510  return pure_saturation.evaluate(iSmolar, _p, _Q, cached_saturation_iL, cached_saturation_iV);
511  }
512  }
513 }
514 CoolPropDbl CoolProp::TabularBackend::calc_umolar(void) {
515  PhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
516  PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
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:
524  throw ValueError("table not selected");
525  }
526  return _HUGE; // not needed, will never be hit, just to make compiler happy
527  } else {
528  if (is_mixture) {
529  // h = u + pv -> u = h - pv
530  CoolPropDbl hmolar = phase_envelope_sat(phase_envelope, iHmolar, iP, _p);
531  CoolPropDbl rhomolar = phase_envelope_sat(phase_envelope, iDmolar, iP, _p);
532  return hmolar - _p / rhomolar;
533  } else {
534  return pure_saturation.evaluate(iUmolar, _p, _Q, cached_saturation_iL, cached_saturation_iV);
535  }
536  }
537 }
538 CoolPropDbl CoolProp::TabularBackend::calc_cpmolar(void) {
539  PhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
540  PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
541  if (using_single_phase_table) {
542  return calc_first_partial_deriv(iHmolar, iT, iP);
543  } else {
544  if (is_mixture) {
545  return phase_envelope_sat(phase_envelope, iCpmolar, iP, _p);
546  } else {
547  return pure_saturation.evaluate(iCpmolar, _p, _Q, cached_saturation_iL, cached_saturation_iV);
548  }
549  }
550 }
551 CoolPropDbl CoolProp::TabularBackend::calc_cvmolar(void) {
552  PhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
553  PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
554  if (using_single_phase_table) {
555  return calc_first_partial_deriv(iUmolar, iT, iDmolar);
556  } else {
557  if (is_mixture) {
558  return phase_envelope_sat(phase_envelope, iCvmolar, iP, _p);
559  } else {
560  return pure_saturation.evaluate(iCvmolar, _p, _Q, cached_saturation_iL, cached_saturation_iV);
561  }
562  }
563 }
564 
565 CoolPropDbl CoolProp::TabularBackend::calc_viscosity(void) {
566  PhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
567  PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
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:
575  throw ValueError("table not selected");
576  }
577  return _HUGE; // not needed, will never be hit, just to make compiler happy
578  } else {
579  if (is_mixture) {
580  return phase_envelope_sat(phase_envelope, iviscosity, iP, _p);
581  } else {
582  return pure_saturation.evaluate(iviscosity, _p, _Q, cached_saturation_iL, cached_saturation_iV);
583  }
584  }
585 }
586 CoolPropDbl CoolProp::TabularBackend::calc_conductivity(void) {
587  PhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
588  PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
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:
596  throw ValueError("table not selected");
597  }
598  return _HUGE; // not needed, will never be hit, just to make compiler happy
599  } else {
600  if (is_mixture) {
601  return phase_envelope_sat(phase_envelope, iconductivity, iP, _p);
602  } else {
603  return pure_saturation.evaluate(iconductivity, _p, _Q, cached_saturation_iL, cached_saturation_iV);
604  }
605  }
606 }
607 CoolPropDbl CoolProp::TabularBackend::calc_speed_sound(void) {
608  PhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
609  PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
610  if (using_single_phase_table) {
611  return sqrt(1 / molar_mass() * cpmolar() / cvmolar() * first_partial_deriv(iP, iDmolar, iT));
612  } else {
613  if (is_mixture) {
614  return phase_envelope_sat(phase_envelope, ispeed_sound, iP, _p);
615  } else {
616  return pure_saturation.evaluate(ispeed_sound, _p, _Q, cached_saturation_iL, cached_saturation_iV);
617  }
618  }
619 }
620 CoolPropDbl CoolProp::TabularBackend::calc_first_partial_deriv(parameters Of, parameters Wrt, parameters Constant) {
621  if (using_single_phase_table) {
622  CoolPropDbl dOf_dx, dOf_dy, dWrt_dx, dWrt_dy, dConstant_dx, dConstant_dy;
623 
624  // If a mass-based parameter is provided, get a conversion factor and change the key to the molar-based key
625  double Of_conversion_factor = 1.0, Wrt_conversion_factor = 1.0, Constant_conversion_factor = 1.0, MM = AS->molar_mass();
626  mass_to_molar(Of, Of_conversion_factor, MM);
627  mass_to_molar(Wrt, Wrt_conversion_factor, MM);
628  mass_to_molar(Constant, Constant_conversion_factor, MM);
629 
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);
638  break;
639  }
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);
647  break;
648  }
649  case SELECTED_NO_TABLE:
650  throw ValueError("table not selected");
651  }
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;
654  } else {
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));
656  }
657 };
658 
659 CoolPropDbl CoolProp::TabularBackend::calc_first_saturation_deriv(parameters Of1, parameters Wrt1) {
660  PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
661  if (AS->get_mole_fractions().size() > 1) {
662  throw ValueError("calc_first_saturation_deriv not available for mixtures");
663  }
664  if (std::abs(_Q) < 1e-6) {
665  return pure_saturation.first_saturation_deriv(Of1, Wrt1, 0, keyed_output(Wrt1), cached_saturation_iL);
666  } else if (std::abs(_Q - 1) < 1e-6) {
667  return pure_saturation.first_saturation_deriv(Of1, Wrt1, 1, keyed_output(Wrt1), cached_saturation_iV);
668  } else {
669  throw ValueError(format("Quality [%Lg] must be either 0 or 1 to within 1 ppm", _Q));
670  }
671 }
672 CoolPropDbl CoolProp::TabularBackend::calc_first_two_phase_deriv(parameters Of, parameters Wrt, parameters Constant) {
673  PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
674  if (Of == iDmolar && Wrt == iHmolar && Constant == iP) {
675  CoolPropDbl rhoL = pure_saturation.evaluate(iDmolar, _p, 0, cached_saturation_iL, cached_saturation_iV);
676  CoolPropDbl rhoV = pure_saturation.evaluate(iDmolar, _p, 1, cached_saturation_iL, cached_saturation_iV);
677  CoolPropDbl hL = pure_saturation.evaluate(iHmolar, _p, 0, cached_saturation_iL, cached_saturation_iV);
678  CoolPropDbl hV = pure_saturation.evaluate(iHmolar, _p, 1, cached_saturation_iL, cached_saturation_iV);
679  return -POW2(rhomolar()) * (1 / rhoV - 1 / rhoL) / (hV - hL);
680  } else if (Of == iDmass && Wrt == iHmass && Constant == iP) {
681  return first_two_phase_deriv(iDmolar, iHmolar, iP) * POW2(molar_mass());
682  } else if (Of == iDmolar && Wrt == iP && Constant == iHmolar) {
683  // v = 1/rho; drhodv = -rho^2; dvdrho = -1/rho^2
684  CoolPropDbl rhoL = pure_saturation.evaluate(iDmolar, _p, 0, cached_saturation_iL, cached_saturation_iV);
685  CoolPropDbl rhoV = pure_saturation.evaluate(iDmolar, _p, 1, cached_saturation_iL, cached_saturation_iV);
686  CoolPropDbl hL = pure_saturation.evaluate(iHmolar, _p, 0, cached_saturation_iL, cached_saturation_iV);
687  CoolPropDbl hV = pure_saturation.evaluate(iHmolar, _p, 1, cached_saturation_iL, cached_saturation_iV);
688  CoolPropDbl dvdrhoL = -1 / POW2(rhoL);
689  CoolPropDbl dvdrhoV = -1 / POW2(rhoV);
690  CoolPropDbl dvL_dp = dvdrhoL * pure_saturation.first_saturation_deriv(iDmolar, iP, 0, _p, cached_saturation_iL);
691  CoolPropDbl dvV_dp = dvdrhoV * pure_saturation.first_saturation_deriv(iDmolar, iP, 1, _p, cached_saturation_iV);
692  CoolPropDbl dhL_dp = pure_saturation.first_saturation_deriv(iHmolar, iP, 0, _p, cached_saturation_iL);
693  CoolPropDbl dhV_dp = pure_saturation.first_saturation_deriv(iHmolar, iP, 1, _p, cached_saturation_iV);
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;
697  } else if (Of == iDmass && Wrt == iP && Constant == iHmass) {
698  return first_two_phase_deriv(iDmolar, iP, iHmolar) * molar_mass();
699  } else {
700  throw ValueError("These inputs are not supported to calc_first_two_phase_deriv");
701  }
702 }
703 
704 CoolPropDbl CoolProp::TabularBackend::calc_first_two_phase_deriv_splined(parameters Of, parameters Wrt, parameters Constant, CoolPropDbl x_end) {
705  // Note: If you need all three values (drho_dh__p, drho_dp__h and rho_spline),
706  // you should calculate drho_dp__h first to avoid duplicate calculations.
707 
708  bool drho_dh__p = false;
709  bool drho_dp__h = false;
710  bool rho_spline = false;
711 
712  if (Of == iDmolar && Wrt == iHmolar && Constant == iP) {
713  drho_dh__p = true;
714  if (_drho_spline_dh__constp) return _drho_spline_dh__constp;
715  } else if (Of == iDmass && Wrt == iHmass && Constant == iP) {
716  return first_two_phase_deriv_splined(iDmolar, iHmolar, iP, x_end) * POW2(molar_mass());
717  } else if (Of == iDmolar && Wrt == iP && Constant == iHmolar) {
718  drho_dp__h = true;
719  if (_drho_spline_dp__consth) return _drho_spline_dp__consth;
720  } else if (Of == iDmass && Wrt == iP && Constant == iHmass) {
721  return first_two_phase_deriv_splined(iDmolar, iP, iHmolar, x_end) * molar_mass();
722  }
723  // Add the special case for the splined density
724  else if (Of == iDmolar && Wrt == iDmolar && Constant == iDmolar) {
725  rho_spline = true;
726  if (_rho_spline) return _rho_spline;
727  } else if (Of == iDmass && Wrt == iDmass && Constant == iDmass) {
728  return first_two_phase_deriv_splined(iDmolar, iDmolar, iDmolar, x_end) * molar_mass();
729  } else {
730  throw ValueError("These inputs are not supported to calc_first_two_phase_deriv");
731  }
732 
733  if (_Q > x_end) {
734  throw ValueError(format("Q [%g] is greater than x_end [%Lg]", _Q, x_end).c_str());
735  }
736  if (_phase != iphase_twophase) {
737  throw ValueError(format("state is not two-phase"));
738  }
739 
740  // TODO: replace AS with a cloned instance of TTSE or BICUBIC, add "clone()" as mandatory function in base class
741  //shared_ptr<TabularBackend>
742  // Liq(this->clone())),
743  // End(this->clone()));
744  //
745  //Liq->specify_phase(iphase_liquid);
746  //Liq->_Q = -1;
747  //Liq->update_DmolarT_direct(SatL->rhomolar(), SatL->T());
748  //End->update(QT_INPUTS, x_end, SatL->T());
749 
750  // Start with quantities needed for all calculations
751  PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
752  CoolPropDbl hL = pure_saturation.evaluate(iHmolar, _p, 0, cached_saturation_iL, cached_saturation_iV);
753  CoolPropDbl hV = pure_saturation.evaluate(iHmolar, _p, 1, cached_saturation_iL, cached_saturation_iV);
754  CoolPropDbl hE = pure_saturation.evaluate(iHmolar, _p, x_end, cached_saturation_iL, cached_saturation_iV);
755 
756  CoolPropDbl dL = pure_saturation.evaluate(iDmolar, _p, 0, cached_saturation_iL, cached_saturation_iV);
757  CoolPropDbl dV = pure_saturation.evaluate(iDmolar, _p, 1, cached_saturation_iL, cached_saturation_iV);
758  CoolPropDbl dE = pure_saturation.evaluate(iDmolar, _p, x_end, cached_saturation_iL, cached_saturation_iV);
759 
760  CoolPropDbl Delta = Q() * (hV - hL);
761  CoolPropDbl Delta_end = hE - hL;
762 
763  CoolPropDbl TL = pure_saturation.evaluate(iT, _p, 0, cached_saturation_iL, cached_saturation_iV);
764  CoolPropDbl TE = pure_saturation.evaluate(iT, _p, x_end, cached_saturation_iL, cached_saturation_iV);
765 
766  // TODO: We cheat here and fake access to a derived class.
767  // Liquid state
768  AS->specify_phase(iphase_liquid);
769  AS->update(DmolarT_INPUTS, dL, TL);
770  CoolPropDbl drho_dh_liq__constp = AS->first_partial_deriv(iDmolar, iHmolar, iP);
771  CoolPropDbl d2rhodhdp_liq = AS->second_partial_deriv(iDmolar, iHmolar, iP, iP, iHmolar);
772  // End of spline
773  AS->specify_phase(iphase_twophase);
774  AS->update(DmolarT_INPUTS, dE, TE);
775  CoolPropDbl drho_dh_end__constp = AS->first_partial_deriv(iDmolar, iHmolar, iP);
776  CoolPropDbl d2rhodhdp_end = AS->second_partial_deriv(iDmolar, iHmolar, iP, iP, iHmolar);
777 
778  // Spline coordinates a, b, c, d
779  CoolPropDbl Abracket = (2 * dL - 2 * dE + Delta_end * (drho_dh_liq__constp + drho_dh_end__constp));
780  CoolPropDbl a = 1 / POW3(Delta_end) * Abracket;
781  CoolPropDbl b = 3 / POW2(Delta_end) * (-dL + dE) - 1 / Delta_end * (drho_dh_end__constp + 2 * drho_dh_liq__constp);
782  CoolPropDbl c = drho_dh_liq__constp;
783  CoolPropDbl d = dL;
784 
785  // Either the spline value or drho/dh|p can be directly evaluated now
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;
790 
791  // It's drho/dp|h
792  // ... calculate some more things
793 
794  // Derivatives *along* the saturation curve using the special internal method
795  CoolPropDbl dhL_dp_sat = pure_saturation.first_saturation_deriv(iHmolar, iP, 0, _p, cached_saturation_iL);
796  CoolPropDbl dhV_dp_sat = pure_saturation.first_saturation_deriv(iHmolar, iP, 1, _p, cached_saturation_iV);
797  CoolPropDbl drhoL_dp_sat = pure_saturation.first_saturation_deriv(iDmolar, iP, 0, _p, cached_saturation_iL);
798  CoolPropDbl drhoV_dp_sat = pure_saturation.first_saturation_deriv(iDmolar, iP, 1, _p, cached_saturation_iV);
799  //CoolPropDbl rhoV = SatV->keyed_output(rho_key);
800  //CoolPropDbl rhoL = SatL->keyed_output(rho_key);
801  CoolPropDbl drho_dp_end = POW2(dE) * (x_end / POW2(dV) * drhoV_dp_sat + (1 - x_end) / POW2(dL) * drhoL_dp_sat);
802 
803  // Faking single-phase
804  //CoolPropDbl drho_dp__consth_liq = Liq->first_partial_deriv(rho_key, p_key, h_key);
805  //CoolPropDbl d2rhodhdp_liq = Liq->second_partial_deriv(rho_key, h_key, p_key, p_key, h_key); // ?
806 
807  // Derivatives at the end point
808  // CoolPropDbl drho_dp__consth_end = End->calc_first_two_phase_deriv(rho_key, p_key, h_key);
809  //CoolPropDbl d2rhodhdp_end = End->calc_second_two_phase_deriv(rho_key, h_key, p_key, p_key, h_key);
810 
811  // Reminder:
812  // Delta = Q()*(hV-hL) = h-hL
813  // Delta_end = x_end*(hV-hL);
814  CoolPropDbl d_Delta_dp__consth = -dhL_dp_sat;
815  CoolPropDbl d_Delta_end_dp__consth = x_end * (dhV_dp_sat - dhL_dp_sat);
816 
817  // First pressure derivative at constant h of the coefficients a,b,c,d
818  // CoolPropDbl Abracket = (2*rho_liq - 2*rho_end + Delta_end * (drho_dh_liq__constp + drho_dh_end));
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);
825  CoolPropDbl dc_dp = d2rhodhdp_liq;
826  CoolPropDbl dd_dp = drhoL_dp_sat;
827 
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;
831 
832  throw ValueError("Something went wrong in TabularBackend::calc_first_two_phase_deriv_splined");
833  return _HUGE;
834 }
835 
836 void CoolProp::TabularBackend::update(CoolProp::input_pairs input_pair, double val1, double val2) {
837 
838  if (get_debug_level() > 0) {
839  std::cout << format("update(%s,%g,%g)\n", get_input_pair_short_desc(input_pair).c_str(), val1, val2);
840  }
841 
842  // Clear cached variables
843  clear();
844 
845  // Convert to mass-based units if necessary
846  CoolPropDbl ld_value1 = val1, ld_value2 = val2;
847  mass_to_molar_inputs(input_pair, ld_value1, ld_value2);
848  val1 = ld_value1;
849  val2 = ld_value2;
850 
851  // Check the tables, build if necessary
852  check_tables();
853 
854  // Flush the cached indices (set to large number)
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();
859 
860  // To start, set quality to value that is impossible
861  _Q = -1000;
862 
863  PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
864  PhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
865  SinglePhaseGriddedTableData& single_phase_logph = dataset->single_phase_logph;
866  SinglePhaseGriddedTableData& single_phase_logpT = dataset->single_phase_logpT;
867 
868  switch (input_pair) {
869  case HmolarP_INPUTS: {
870  _hmolar = val1;
871  _p = val2;
872  if (!single_phase_logph.native_inputs_are_in_range(_hmolar, _p)) {
873  // Use the AbstractState instance
874  using_single_phase_table = false;
875  if (get_debug_level() > 5) {
876  std::cout << "inputs are not in range";
877  }
878  throw ValueError(format("inputs are not in range, hmolar=%Lg, p=%Lg", static_cast<CoolPropDbl>(_hmolar), _p));
879  } else {
880  using_single_phase_table = true; // Use the table !
881  std::size_t iL = std::numeric_limits<std::size_t>::max(), iV = std::numeric_limits<std::size_t>::max(), iclosest = 0;
882  CoolPropDbl hL = 0, hV = 0;
883  SimpleState closest_state;
884  bool is_two_phase = false;
885  // If phase is imposed, use it, but only if it's single phase:
886  // - Imposed two phase still needs to determine saturation limits by calling pure_saturation.is_inside().
887  // - There's no speed increase to be gained by imposing two phase.
888  if ((imposed_phase_index == iphase_not_imposed) || (imposed_phase_index == iphase_twophase)) {
889  if (is_mixture) {
890  is_two_phase = PhaseEnvelopeRoutines::is_inside(phase_envelope, iP, _p, iHmolar, _hmolar, iclosest, closest_state);
891  } else {
892  is_two_phase = pure_saturation.is_inside(iP, _p, iHmolar, _hmolar, iL, iV, hL, hV);
893  }
894  }
895  // Phase determined or imposed, now interpolate results
896  if (is_two_phase) {
897  using_single_phase_table = false;
898  _Q = (static_cast<double>(_hmolar) - hL) / (hV - hL);
899  if (!is_in_closed_range(0.0, 1.0, static_cast<double>(_Q))) {
900  throw ValueError(
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));
902  } else {
903  cached_saturation_iL = iL;
904  cached_saturation_iV = iV;
905  _phase = iphase_twophase;
906  }
907  } else {
908  selected_table = SELECTED_PH_TABLE;
909  // Find and cache the indices i, j
910  find_native_nearest_good_indices(single_phase_logph, dataset->coeffs_ph, _hmolar, _p, cached_single_phase_i,
911  cached_single_phase_j);
912  // Recalculate the phase
913  recalculate_singlephase_phase();
914  }
915  }
916  break;
917  }
918  case PT_INPUTS: {
919  _p = val1;
920  _T = val2;
921  if (!single_phase_logpT.native_inputs_are_in_range(_T, _p)) {
922  // Use the AbstractState instance
923  using_single_phase_table = false;
924  if (get_debug_level() > 5) {
925  std::cout << "inputs are not in range";
926  }
927  throw ValueError(format("inputs are not in range, p=%g Pa, T=%g K", _p, _T));
928  } else {
929  using_single_phase_table = true; // Use the table !
930  std::size_t iL = 0, iV = 0, iclosest = 0;
931  CoolPropDbl TL = 0, TV = 0;
932  SimpleState closest_state;
933  bool is_two_phase = false;
934  // Phase is imposed, use it
935  if (imposed_phase_index != iphase_not_imposed) {
936  is_two_phase = (imposed_phase_index == iphase_twophase);
937  } else {
938  if (is_mixture) {
939  is_two_phase = PhaseEnvelopeRoutines::is_inside(phase_envelope, iP, _p, iT, _T, iclosest, closest_state);
940  } else {
941  is_two_phase = pure_saturation.is_inside(iP, _p, iT, _T, iL, iV, TL, TV);
942  }
943  }
944  if (is_two_phase) {
945  using_single_phase_table = false;
946  throw ValueError(format("P,T with TTSE cannot be two-phase for now"));
947  } else {
948  selected_table = SELECTED_PT_TABLE;
949  // Find and cache the indices i, j
950  find_native_nearest_good_indices(single_phase_logpT, dataset->coeffs_pT, _T, _p, cached_single_phase_i, cached_single_phase_j);
951 
952  if (imposed_phase_index != iphase_not_imposed) {
953  double rhoc = rhomolar_critical();
954  if (imposed_phase_index == iphase_liquid && cached_single_phase_i > 0) {
955  // We want a liquid solution, but we got a vapor solution
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) {
959  // Bump to lower temperature
960  cached_single_phase_i--;
961  }
962  double rho = evaluate_single_phase_pT(iDmolar, cached_single_phase_i, cached_single_phase_j);
963  if (rho < rhoc) {
964  // Didn't work
965  throw ValueError("Bump unsuccessful");
966  } else {
967  _rhomolar = rho;
968  }
969  }
970  } else if ((imposed_phase_index == iphase_gas || imposed_phase_index == iphase_supercritical_gas)
971  && cached_single_phase_i > 0) {
972 
973  // We want a gas solution, but we got a liquid solution
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) {
977  // Bump to lower temperature
978  cached_single_phase_i++;
979  }
980  double rho = evaluate_single_phase_pT(iDmolar, cached_single_phase_i, cached_single_phase_j);
981  if (rho > rhoc) {
982  // Didn't work
983  throw ValueError("Bump unsuccessful");
984  } else {
985  _rhomolar = rho;
986  }
987  }
988  }
989  } else {
990 
991  // If p < pc, you might be getting a liquid solution when you want a vapor solution or vice versa
992  // if you are very close to the saturation curve, so we figure out what the saturation temperature
993  // is for the given pressure
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) {
999  if (_T < Ts) {
1000  if (cached_single_phase_i == 0) {
1001  throw ValueError(format("P, T are near saturation, but cannot move the cell to the left"));
1002  }
1003  // It's liquid, move the cell to the left
1004  cached_single_phase_i--;
1005  } else {
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"));
1008  }
1009  // It's vapor, move to the right
1010  cached_single_phase_i++;
1011  }
1012  }
1013  }
1014  }
1015  // Recalculate the phase
1016  recalculate_singlephase_phase();
1017  }
1018  }
1019  break;
1020  }
1021  case PUmolar_INPUTS:
1022  case PSmolar_INPUTS:
1023  case DmolarP_INPUTS: {
1024  CoolPropDbl otherval;
1025  parameters otherkey;
1026  switch (input_pair) {
1027  case PUmolar_INPUTS:
1028  _p = val1;
1029  _umolar = val2;
1030  otherval = val2;
1031  otherkey = iUmolar;
1032  break;
1033  case PSmolar_INPUTS:
1034  _p = val1;
1035  _smolar = val2;
1036  otherval = val2;
1037  otherkey = iSmolar;
1038  break;
1039  case DmolarP_INPUTS:
1040  _rhomolar = val1;
1041  _p = val2;
1042  otherval = val1;
1043  otherkey = iDmolar;
1044  break;
1045  default:
1046  throw ValueError("Bad (impossible) pair");
1047  }
1048 
1049  using_single_phase_table = true; // Use the table (or first guess is that it is single-phase)!
1050  std::size_t iL = std::numeric_limits<std::size_t>::max(), iV = std::numeric_limits<std::size_t>::max();
1051  CoolPropDbl zL = 0, zV = 0;
1052  std::size_t iclosest = 0;
1053  SimpleState closest_state;
1054  bool is_two_phase = false;
1055  // If phase is imposed, use it, but only if it's single phase:
1056  // - Imposed two phase still needs to determine saturation limits by calling is_inside().
1057  // - There's no speed increase to be gained by imposing two phase.
1058  if ((imposed_phase_index == iphase_not_imposed) || (imposed_phase_index == iphase_twophase)) {
1059  if (is_mixture) {
1060  is_two_phase = PhaseEnvelopeRoutines::is_inside(phase_envelope, iP, _p, otherkey, otherval, iclosest, closest_state);
1061  if (is_two_phase) {
1062  std::vector<std::pair<std::size_t, std::size_t>> intersect =
1063  PhaseEnvelopeRoutines::find_intersections(phase_envelope, iP, _p);
1064  if (intersect.size() < 2) {
1065  throw ValueError(format("p [%g Pa] is not within phase envelope", _p));
1066  }
1067  iV = intersect[0].first;
1068  iL = intersect[1].first;
1069  zL = PhaseEnvelopeRoutines::evaluate(phase_envelope, otherkey, iP, _p, iL);
1070  zV = PhaseEnvelopeRoutines::evaluate(phase_envelope, otherkey, iP, _p, iV);
1071  }
1072  } else {
1073  is_two_phase = pure_saturation.is_inside(iP, _p, otherkey, otherval, iL, iV, zL, zV);
1074  }
1075  }
1076  // Phase determined or imposed, now interpolate results
1077  if (is_two_phase) {
1078  using_single_phase_table = false;
1079  if (otherkey == iDmolar) {
1080  _Q = (1 / otherval - 1 / zL) / (1 / zV - 1 / zL);
1081  } else {
1082  _Q = (otherval - zL) / (zV - zL);
1083  }
1084  if (!is_in_closed_range(0.0, 1.0, static_cast<double>(_Q))) {
1085 
1086  throw ValueError(format("vapor quality is not in (0,1) for %s: %g p: %g", get_parameter_information(otherkey, "short").c_str(),
1087  otherval, static_cast<double>(_p)));
1088  } else if (!is_mixture) {
1089  cached_saturation_iL = iL;
1090  cached_saturation_iV = iV;
1091  }
1092  _phase = iphase_twophase;
1093  } else {
1094  selected_table = SELECTED_PH_TABLE;
1095  // Find and cache the indices i, j
1096  find_nearest_neighbor(single_phase_logph, dataset->coeffs_ph, iP, _p, otherkey, otherval, cached_single_phase_i,
1097  cached_single_phase_j);
1098  // Now find hmolar given P, X for X in Smolar, Umolar, Dmolar
1099  invert_single_phase_x(single_phase_logph, dataset->coeffs_ph, otherkey, otherval, _p, cached_single_phase_i, cached_single_phase_j);
1100  // Recalculate the phase
1101  recalculate_singlephase_phase();
1102  }
1103  break;
1104  }
1105  case SmolarT_INPUTS:
1106  case DmolarT_INPUTS: {
1107  CoolPropDbl otherval;
1108  parameters otherkey;
1109  switch (input_pair) {
1110  case SmolarT_INPUTS:
1111  _smolar = val1;
1112  _T = val2;
1113  otherval = val1;
1114  otherkey = iSmolar;
1115  break;
1116  case DmolarT_INPUTS:
1117  _rhomolar = val1;
1118  _T = val2;
1119  otherval = val1;
1120  otherkey = iDmolar;
1121  break;
1122  default:
1123  throw ValueError("Bad (impossible) pair");
1124  }
1125 
1126  using_single_phase_table = true; // Use the table (or first guess is that it is single-phase)!
1127  std::size_t iL = std::numeric_limits<std::size_t>::max(), iV = std::numeric_limits<std::size_t>::max();
1128  CoolPropDbl zL = 0, zV = 0;
1129  std::size_t iclosest = 0;
1130  SimpleState closest_state;
1131  bool is_two_phase = false;
1132  // If phase is imposed, use it, but only if it's single phase:
1133  // - Imposed two phase still needs to determine saturation limits by calling is_inside().
1134  // - There's no speed increase to be gained by imposing two phase.
1135  if ((imposed_phase_index == iphase_not_imposed) || (imposed_phase_index == iphase_twophase)) {
1136  if (is_mixture) {
1137  is_two_phase = PhaseEnvelopeRoutines::is_inside(phase_envelope, iT, _T, otherkey, otherval, iclosest, closest_state);
1138  if (is_two_phase) {
1139  std::vector<std::pair<std::size_t, std::size_t>> intersect =
1140  PhaseEnvelopeRoutines::find_intersections(phase_envelope, iT, _T);
1141  if (intersect.size() < 2) {
1142  throw ValueError(format("T [%g K] is not within phase envelope", _T));
1143  }
1144  iV = intersect[0].first;
1145  iL = intersect[1].first;
1146  zL = PhaseEnvelopeRoutines::evaluate(phase_envelope, otherkey, iT, _T, iL);
1147  zV = PhaseEnvelopeRoutines::evaluate(phase_envelope, otherkey, iT, _T, iV);
1148  }
1149  } else {
1150  is_two_phase = pure_saturation.is_inside(iT, _T, otherkey, otherval, iL, iV, zL, zV);
1151  }
1152  }
1153  if (is_two_phase) {
1154  using_single_phase_table = false;
1155  if (otherkey == iDmolar) {
1156  _Q = (1 / otherval - 1 / zL) / (1 / zV - 1 / zL);
1157  } else {
1158  _Q = (otherval - zL) / (zV - zL);
1159  }
1160  if (!is_in_closed_range(0.0, 1.0, static_cast<double>(_Q))) {
1161  throw ValueError(format("vapor quality is not in (0,1) for %s: %g T: %g", get_parameter_information(otherkey, "short").c_str(),
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);
1167  } else {
1168  // Mixture
1169  std::vector<std::pair<std::size_t, std::size_t>> intersect = PhaseEnvelopeRoutines::find_intersections(phase_envelope, iT, _T);
1170  if (intersect.size() < 2) {
1171  throw ValueError(format("T [%g K] is not within phase envelope", _T));
1172  }
1173  iV = intersect[0].first;
1174  iL = intersect[1].first;
1175  CoolPropDbl pL = PhaseEnvelopeRoutines::evaluate(phase_envelope, iP, iT, _T, iL);
1176  CoolPropDbl pV = PhaseEnvelopeRoutines::evaluate(phase_envelope, iP, iT, _T, iV);
1177  _p = _Q * pV + (1 - _Q) * pL;
1178  }
1179  } else {
1180  selected_table = SELECTED_PT_TABLE;
1181  // Find and cache the indices i, j
1182  find_nearest_neighbor(single_phase_logpT, dataset->coeffs_pT, iT, _T, otherkey, otherval, cached_single_phase_i,
1183  cached_single_phase_j);
1184  // Now find the y variable (Dmolar or Smolar in this case)
1185  invert_single_phase_y(single_phase_logpT, dataset->coeffs_pT, otherkey, otherval, _T, cached_single_phase_i, cached_single_phase_j);
1186  // Recalculate the phase
1187  recalculate_singlephase_phase();
1188  }
1189  break;
1190  }
1191  case PQ_INPUTS: {
1192  std::size_t iL = 0, iV = 0;
1193  _p = val1;
1194  _Q = val2;
1195  using_single_phase_table = false;
1196  if (!is_in_closed_range(0.0, 1.0, static_cast<double>(_Q))) {
1197  throw ValueError(format("vapor quality [%g] is not in (0,1)", static_cast<double>(val2)));
1198  } else {
1199  CoolPropDbl TL = _HUGE, TV = _HUGE;
1200  if (is_mixture) {
1201  std::vector<std::pair<std::size_t, std::size_t>> intersect = PhaseEnvelopeRoutines::find_intersections(phase_envelope, iP, _p);
1202  if (intersect.size() < 2) {
1203  throw ValueError(format("p [%g Pa] is not within phase envelope", _p));
1204  }
1205  iV = intersect[0].first;
1206  iL = intersect[1].first;
1207  TL = PhaseEnvelopeRoutines::evaluate(phase_envelope, iT, iP, _p, iL);
1208  TV = PhaseEnvelopeRoutines::evaluate(phase_envelope, iT, iP, _p, iV);
1209  } else {
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");
1213  }
1214  }
1215  _T = _Q * TV + (1 - _Q) * TL;
1216  cached_saturation_iL = iL;
1217  cached_saturation_iV = iV;
1218  _phase = iphase_twophase;
1219  }
1220  break;
1221  }
1222  case QT_INPUTS: {
1223  std::size_t iL = 0, iV = 0;
1224  _Q = val1;
1225  _T = val2;
1226 
1227  using_single_phase_table = false;
1228  if (!is_in_closed_range(0.0, 1.0, static_cast<double>(_Q))) {
1229  throw ValueError(format("vapor quality [%g] is not in (0,1)", static_cast<double>(val1)));
1230  } else {
1231  CoolPropDbl pL = _HUGE, pV = _HUGE;
1232  if (is_mixture) {
1233  std::vector<std::pair<std::size_t, std::size_t>> intersect = PhaseEnvelopeRoutines::find_intersections(phase_envelope, iT, _T);
1234  if (intersect.size() < 2) {
1235  throw ValueError(format("T [%g K] is not within phase envelope", _T));
1236  }
1237  iV = intersect[0].first;
1238  iL = intersect[1].first;
1239  pL = PhaseEnvelopeRoutines::evaluate(phase_envelope, iP, iT, _T, iL);
1240  pV = PhaseEnvelopeRoutines::evaluate(phase_envelope, iP, iT, _T, iV);
1241  } else {
1242  pure_saturation.is_inside(iT, _T, iQ, _Q, iL, iV, pL, pV);
1243  }
1244  _p = _Q * pV + (1 - _Q) * pL;
1245  cached_saturation_iL = iL;
1246  cached_saturation_iV = iV;
1247  _phase = iphase_twophase;
1248  }
1249  break;
1250  }
1251  default:
1252  throw ValueError("Sorry, but this set of inputs is not supported for Tabular backend");
1253  }
1254 }
1255 
1256 void CoolProp::TabularDataSet::write_tables(const std::string& path_to_tables) {
1257  make_dirs(path_to_tables);
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");
1262 }
1263 
1264 void CoolProp::TabularDataSet::load_tables(const std::string& path_to_tables, shared_ptr<CoolProp::AbstractState>& AS) {
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;
1275  if (get_debug_level() > 0) {
1276  std::cout << "Tables loaded" << std::endl;
1277  }
1278 };
1279 
1280 void CoolProp::TabularDataSet::build_tables(shared_ptr<CoolProp::AbstractState>& AS) {
1281  // Pure or pseudo-pure fluid
1282  if (AS->get_mole_fractions().size() == 1) {
1283  pure_saturation.build(AS);
1284  } else {
1285  // Call function to actually construct the phase envelope
1286  AS->build_phase_envelope("");
1287  // Copy constructed phase envelope into this class
1288  PhaseEnvelopeData PED = AS->get_phase_envelope_data();
1289  // Convert into packable form
1290  phase_envelope.copy_from_nonpackable(PED);
1291  // Resize so that it will load properly
1292  pure_saturation.resize(pure_saturation.N);
1293  }
1294  single_phase_logph.build(AS);
1295  single_phase_logpT.build(AS);
1296  tables_loaded = true;
1297 }
1298 
1300 CoolProp::TabularDataSet* CoolProp::TabularDataLibrary::get_set_of_tables(shared_ptr<AbstractState>& AS, bool& loaded) {
1301  const std::string path = path_to_tables(AS);
1302  // Try to find tabular set if it is already loaded
1303  std::map<std::string, TabularDataSet>::iterator it = data.find(path);
1304  // It is already in the map, return it
1305  if (it != data.end()) {
1306  loaded = it->second.tables_loaded;
1307  return &(it->second);
1308  }
1309  // It is not in the map, build it
1310  else {
1311  TabularDataSet set;
1312  data.insert(std::pair<std::string, TabularDataSet>(path, set));
1313  TabularDataSet& dataset = data[path];
1314  try {
1315  if (!dataset.tables_loaded) {
1316  dataset.load_tables(path, AS);
1317  }
1318  loaded = true;
1319  } catch (std::exception&) {
1320  loaded = false;
1321  }
1322  return &(dataset);
1323  }
1324 }
1325 
1326 void CoolProp::TabularDataSet::build_coeffs(SinglePhaseGriddedTableData& table, std::vector<std::vector<CellCoeffs>>& coeffs) {
1327  if (!coeffs.empty()) {
1328  return;
1329  }
1330  const bool debug = get_debug_level() > 5 || false;
1331  const int param_count = 6;
1332  parameters param_list[param_count] = {iDmolar, iT, iSmolar, iHmolar, iP, iUmolar};
1333  std::vector<std::vector<double>>*f = NULL, *fx = NULL, *fy = NULL, *fxy = NULL;
1334 
1335  clock_t t1 = clock();
1336 
1337  // Resize the coefficient structures
1338  coeffs.resize(table.Nx - 1, std::vector<CellCoeffs>(table.Ny - 1));
1339 
1340  int valid_cell_count = 0;
1341  for (std::size_t k = 0; k < param_count; ++k) {
1342  parameters param = param_list[k];
1343  if (param == table.xkey || param == table.ykey) {
1344  continue;
1345  } // Skip tables that match either of the input variables
1346 
1347  switch (param) {
1348  case iT:
1349  f = &(table.T);
1350  fx = &(table.dTdx);
1351  fy = &(table.dTdy);
1352  fxy = &(table.d2Tdxdy);
1353  break;
1354  case iP:
1355  f = &(table.p);
1356  fx = &(table.dpdx);
1357  fy = &(table.dpdy);
1358  fxy = &(table.d2pdxdy);
1359  break;
1360  case iDmolar:
1361  f = &(table.rhomolar);
1362  fx = &(table.drhomolardx);
1363  fy = &(table.drhomolardy);
1364  fxy = &(table.d2rhomolardxdy);
1365  break;
1366  case iSmolar:
1367  f = &(table.smolar);
1368  fx = &(table.dsmolardx);
1369  fy = &(table.dsmolardy);
1370  fxy = &(table.d2smolardxdy);
1371  break;
1372  case iHmolar:
1373  f = &(table.hmolar);
1374  fx = &(table.dhmolardx);
1375  fy = &(table.dhmolardy);
1376  fxy = &(table.d2hmolardxdy);
1377  break;
1378  case iUmolar:
1379  f = &(table.umolar);
1380  fx = &(table.dumolardx);
1381  fy = &(table.dumolardy);
1382  fxy = &(table.d2umolardxdy);
1383  break;
1384  default:
1385  throw ValueError("Invalid variable type to build_coeffs");
1386  }
1387  for (std::size_t i = 0; i < table.Nx - 1; ++i) // -1 since we have one fewer cells than nodes
1388  {
1389  for (std::size_t j = 0; j < table.Ny - 1; ++j) // -1 since we have one fewer cells than nodes
1390  {
1391  if (ValidNumber((*f)[i][j]) && ValidNumber((*f)[i + 1][j]) && ValidNumber((*f)[i][j + 1]) && ValidNumber((*f)[i + 1][j + 1])) {
1392 
1393  // This will hold the scaled f values for the cell
1394  Eigen::Matrix<double, 16, 1> F;
1395  // The output values (do not require scaling
1396  F(0) = (*f)[i][j];
1397  F(1) = (*f)[i + 1][j];
1398  F(2) = (*f)[i][j + 1];
1399  F(3) = (*f)[i + 1][j + 1];
1400  // Scaling parameter
1401  // d(f)/dxhat = df/dx * dx/dxhat, where xhat = (x-x_i)/(x_{i+1}-x_i)
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;
1408  // Scaling parameter
1409  // d(f)/dyhat = df/dy * dy/dyhat, where yhat = (y-y_j)/(y_{j+1}-y_j)
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;
1416  // Cross derivatives are doubly scaled following the examples above
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;
1421  // Calculate the alpha coefficients
1422  Eigen::MatrixXd alpha = Ainv.transpose() * F; // 16x1; Watch out for the transpose!
1423  std::vector<double> valpha = eigen_to_vec1D(alpha);
1424  coeffs[i][j].set(param, valpha);
1425  coeffs[i][j].set_valid();
1426  valid_cell_count++;
1427  } else {
1428  coeffs[i][j].set_invalid();
1429  }
1430  }
1431  }
1432  double elapsed = (clock() - t1) / ((double)CLOCKS_PER_SEC);
1433  if (debug) {
1434  std::cout << format("Calculated bicubic coefficients for %d good cells in %g sec.\n", valid_cell_count, elapsed);
1435  }
1436  std::size_t remap_count = 0;
1437  // Now find invalid cells and give them pointers to a neighboring cell that works
1438  for (std::size_t i = 0; i < table.Nx - 1; ++i) // -1 since we have one fewer cells than nodes
1439  {
1440  for (std::size_t j = 0; j < table.Ny - 1; ++j) // -1 since we have one fewer cells than nodes
1441  {
1442  // Not a valid cell
1443  if (!coeffs[i][j].valid()) {
1444  // Offsets that we are going to try in order (left, right, top, bottom, diagonals)
1445  int xoffsets[] = {-1, 1, 0, 0, -1, 1, 1, -1};
1446  int yoffsets[] = {0, 0, 1, -1, -1, -1, 1, 1};
1447  // Length of offset
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);
1454  remap_count++;
1455  if (debug) {
1456  std::cout << format("Mapping %d,%d to %d,%d\n", i, j, iplus, jplus);
1457  }
1458  break;
1459  }
1460  }
1461  }
1462  }
1463  }
1464  if (debug) {
1465  std::cout << format("Remapped %d cells\n", remap_count);
1466  }
1467  }
1468 }
1469 
1470 # if defined(ENABLE_CATCH)
1471 # include <catch2/catch_all.hpp>
1472 
1473 // Defined global so we only load once
1474 static shared_ptr<CoolProp::AbstractState> ASHEOS, ASTTSE, ASBICUBIC;
1475 
1476 /* Use a fixture so that loading of the tables from memory only happens once in the initializer */
1477 class TabularFixture
1478 {
1479  public:
1480  TabularFixture() {}
1481  void setup() {
1482  if (ASHEOS.get() == NULL) {
1483  ASHEOS.reset(CoolProp::AbstractState::factory("HEOS", "Water"));
1484  }
1485  if (ASTTSE.get() == NULL) {
1486  ASTTSE.reset(CoolProp::AbstractState::factory("TTSE&HEOS", "Water"));
1487  }
1488  if (ASBICUBIC.get() == NULL) {
1489  ASBICUBIC.reset(CoolProp::AbstractState::factory("BICUBIC&HEOS", "Water"));
1490  }
1491  }
1492 };
1493 TEST_CASE_METHOD(TabularFixture, "Tests for tabular backends with water", "[Tabular]") {
1494  SECTION("first_saturation_deriv invalid quality") {
1495  setup();
1496  ASBICUBIC->update(CoolProp::PQ_INPUTS, 101325, 0.5);
1497  CHECK_THROWS(ASBICUBIC->first_saturation_deriv(CoolProp::iP, CoolProp::iT));
1498  ASTTSE->update(CoolProp::PQ_INPUTS, 101325, 0.5);
1499  CHECK_THROWS(ASTTSE->first_saturation_deriv(CoolProp::iP, CoolProp::iT));
1500  }
1501 
1502  SECTION("first_saturation_deriv dp/dT") {
1503  setup();
1504  ASHEOS->update(CoolProp::PQ_INPUTS, 101325, 1);
1505  CoolPropDbl expected = ASHEOS->first_saturation_deriv(CoolProp::iP, CoolProp::iT);
1506  ASTTSE->update(CoolProp::PQ_INPUTS, 101325, 1);
1507  CoolPropDbl actual_TTSE = ASTTSE->first_saturation_deriv(CoolProp::iP, CoolProp::iT);
1508  ASTTSE->update(CoolProp::PQ_INPUTS, 101325, 1);
1509  CoolPropDbl actual_BICUBIC = ASTTSE->first_saturation_deriv(CoolProp::iP, CoolProp::iT);
1510  CAPTURE(expected);
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);
1515  }
1516  SECTION("first_saturation_deriv dDmolar/dP") {
1517  setup();
1518  ASHEOS->update(CoolProp::PQ_INPUTS, 101325, 1);
1519  CoolPropDbl expected = ASHEOS->first_saturation_deriv(CoolProp::iDmolar, CoolProp::iP);
1520  ASTTSE->update(CoolProp::PQ_INPUTS, 101325, 1);
1521  CoolPropDbl actual_TTSE = ASTTSE->first_saturation_deriv(CoolProp::iDmolar, CoolProp::iP);
1522  ASBICUBIC->update(CoolProp::PQ_INPUTS, 101325, 1);
1523  CoolPropDbl actual_BICUBIC = ASBICUBIC->first_saturation_deriv(CoolProp::iDmolar, CoolProp::iP);
1524  CAPTURE(expected);
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);
1529  }
1530  SECTION("first_saturation_deriv dDmass/dP") {
1531  setup();
1532  ASHEOS->update(CoolProp::PQ_INPUTS, 101325, 1);
1533  CoolPropDbl expected = ASHEOS->first_saturation_deriv(CoolProp::iDmass, CoolProp::iP);
1534  ASTTSE->update(CoolProp::PQ_INPUTS, 101325, 1);
1535  CoolPropDbl actual_TTSE = ASTTSE->first_saturation_deriv(CoolProp::iDmass, CoolProp::iP);
1536  ASBICUBIC->update(CoolProp::PQ_INPUTS, 101325, 1);
1537  CoolPropDbl actual_BICUBIC = ASBICUBIC->first_saturation_deriv(CoolProp::iDmass, CoolProp::iP);
1538  CAPTURE(expected);
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);
1543  }
1544  SECTION("first_saturation_deriv dHmass/dP") {
1545  setup();
1546  ASHEOS->update(CoolProp::PQ_INPUTS, 101325, 1);
1547  CoolPropDbl expected = ASHEOS->first_saturation_deriv(CoolProp::iHmass, CoolProp::iP);
1548  ASTTSE->update(CoolProp::PQ_INPUTS, 101325, 1);
1549  CoolPropDbl actual_TTSE = ASTTSE->first_saturation_deriv(CoolProp::iHmass, CoolProp::iP);
1550  ASBICUBIC->update(CoolProp::PQ_INPUTS, 101325, 1);
1551  CoolPropDbl actual_BICUBIC = ASBICUBIC->first_saturation_deriv(CoolProp::iHmass, CoolProp::iP);
1552  CAPTURE(expected);
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);
1557  }
1558  SECTION("first_saturation_deriv dHmass/dP w/ QT as inputs") {
1559  setup();
1560  ASHEOS->update(CoolProp::QT_INPUTS, 1, 370);
1561  CoolPropDbl expected = ASHEOS->first_saturation_deriv(CoolProp::iHmass, CoolProp::iP);
1562  ASTTSE->update(CoolProp::QT_INPUTS, 1, 370);
1563  CoolPropDbl actual_TTSE = ASTTSE->first_saturation_deriv(CoolProp::iHmass, CoolProp::iP);
1564  ASBICUBIC->update(CoolProp::QT_INPUTS, 1, 370);
1565  CoolPropDbl actual_BICUBIC = ASBICUBIC->first_saturation_deriv(CoolProp::iHmass, CoolProp::iP);
1566  CAPTURE(expected);
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);
1571  }
1572  SECTION("first_two_phase_deriv dDmolar/dP|Hmolar") {
1573  setup();
1574  ASHEOS->update(CoolProp::PQ_INPUTS, 101325, 0.1);
1575  CoolPropDbl expected = ASHEOS->first_two_phase_deriv(CoolProp::iDmolar, CoolProp::iP, CoolProp::iHmolar);
1576  ASTTSE->update(CoolProp::PQ_INPUTS, 101325, 0.1);
1577  CoolPropDbl actual_TTSE = ASTTSE->first_two_phase_deriv(CoolProp::iDmolar, CoolProp::iP, CoolProp::iHmolar);
1578  ASBICUBIC->update(CoolProp::PQ_INPUTS, 101325, 0.1);
1579  CoolPropDbl actual_BICUBIC = ASBICUBIC->first_two_phase_deriv(CoolProp::iDmolar, CoolProp::iP, CoolProp::iHmolar);
1580  CAPTURE(expected);
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);
1585  }
1586  SECTION("first_two_phase_deriv dDmass/dP|Hmass") {
1587  setup();
1588  ASHEOS->update(CoolProp::PQ_INPUTS, 101325, 0.1);
1589  CoolPropDbl expected = ASHEOS->first_two_phase_deriv(CoolProp::iDmass, CoolProp::iP, CoolProp::iHmass);
1590  ASTTSE->update(CoolProp::PQ_INPUTS, 101325, 0.1);
1591  CoolPropDbl actual_TTSE = ASTTSE->first_two_phase_deriv(CoolProp::iDmass, CoolProp::iP, CoolProp::iHmass);
1592  ASBICUBIC->update(CoolProp::PQ_INPUTS, 101325, 0.1);
1593  CoolPropDbl actual_BICUBIC = ASBICUBIC->first_two_phase_deriv(CoolProp::iDmass, CoolProp::iP, CoolProp::iHmass);
1594  CAPTURE(expected);
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);
1599  }
1600  SECTION("first_two_phase_deriv dDmass/dHmass|P") {
1601  setup();
1602  ASHEOS->update(CoolProp::PQ_INPUTS, 101325, 0.1);
1603  CoolPropDbl expected = ASHEOS->first_two_phase_deriv(CoolProp::iDmass, CoolProp::iHmass, CoolProp::iP);
1604  ASTTSE->update(CoolProp::PQ_INPUTS, 101325, 0.1);
1605  CoolPropDbl actual_TTSE = ASTTSE->first_two_phase_deriv(CoolProp::iDmass, CoolProp::iHmass, CoolProp::iP);
1606  ASBICUBIC->update(CoolProp::PQ_INPUTS, 101325, 0.1);
1607  CoolPropDbl actual_BICUBIC = ASBICUBIC->first_two_phase_deriv(CoolProp::iDmass, CoolProp::iHmass, CoolProp::iP);
1608  CAPTURE(expected);
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);
1613  }
1614  SECTION("first_partial_deriv dHmass/dT|P") {
1615  setup();
1616  ASHEOS->update(CoolProp::PT_INPUTS, 101325, 300);
1617  double expected = ASHEOS->cpmass();
1618  ASTTSE->update(CoolProp::PT_INPUTS, 101325, 300);
1619  double dhdT_TTSE = ASTTSE->first_partial_deriv(CoolProp::iHmass, CoolProp::iT, CoolProp::iP);
1620  ASBICUBIC->update(CoolProp::PT_INPUTS, 101325, 300);
1621  double dhdT_BICUBIC = ASBICUBIC->first_partial_deriv(CoolProp::iHmass, CoolProp::iT, CoolProp::iP);
1622  CAPTURE(expected);
1623  CAPTURE(dhdT_TTSE);
1624  CAPTURE(dhdT_BICUBIC);
1625  CHECK(std::abs((expected - dhdT_TTSE) / expected) < 1e-4);
1626  CHECK(std::abs((expected - dhdT_BICUBIC) / expected) < 1e-4);
1627  }
1628  SECTION("first_partial_deriv dHmolar/dT|P") {
1629  setup();
1630  ASHEOS->update(CoolProp::PT_INPUTS, 101325, 300);
1631  double expected = ASHEOS->cpmolar();
1632  ASTTSE->update(CoolProp::PT_INPUTS, 101325, 300);
1633  double dhdT_TTSE = ASTTSE->first_partial_deriv(CoolProp::iHmolar, CoolProp::iT, CoolProp::iP);
1634  ASBICUBIC->update(CoolProp::PT_INPUTS, 101325, 300);
1635  double dhdT_BICUBIC = ASBICUBIC->first_partial_deriv(CoolProp::iHmolar, CoolProp::iT, CoolProp::iP);
1636  CAPTURE(expected);
1637  CAPTURE(dhdT_TTSE);
1638  CAPTURE(dhdT_BICUBIC);
1639  CHECK(std::abs((expected - dhdT_TTSE) / expected) < 1e-4);
1640  CHECK(std::abs((expected - dhdT_BICUBIC) / expected) < 1e-4);
1641  }
1642  SECTION("check isentropic process") {
1643  setup();
1644  double T0 = 300;
1645  double p0 = 1e5;
1646  double p1 = 1e6;
1647  ASHEOS->update(CoolProp::PT_INPUTS, p0, 300);
1648  double s0 = ASHEOS->smolar();
1649  ASHEOS->update(CoolProp::PSmolar_INPUTS, p1, s0);
1650  double expected = ASHEOS->T();
1651  ASTTSE->update(CoolProp::PSmolar_INPUTS, p1, s0);
1652  double actual_TTSE = ASTTSE->T();
1653  ASBICUBIC->update(CoolProp::PSmolar_INPUTS, p1, s0);
1654  double actual_BICUBIC = ASBICUBIC->T();
1655  CAPTURE(expected);
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);
1660  }
1661  SECTION("check D=1 mol/m3, T=500 K inputs") {
1662  setup();
1663  double d = 1;
1664  CAPTURE(d);
1665  ASHEOS->update(CoolProp::DmolarT_INPUTS, d, 500);
1666  double expected = ASHEOS->p();
1667  ASTTSE->update(CoolProp::DmolarT_INPUTS, d, 500);
1668  double actual_TTSE = ASTTSE->p();
1669  ASBICUBIC->update(CoolProp::DmolarT_INPUTS, d, 500);
1670  double actual_BICUBIC = ASBICUBIC->p();
1671  CAPTURE(expected);
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);
1676  }
1677 }
1678 # endif // ENABLE_CATCH
1679 
1680 #endif // !defined(NO_TABULAR_BACKENDS)
Generated on Mon Jan 15 2024 04:22:06 for CoolProp by doxygen 1.9.1