CoolProp  6.6.0
An open-source fluid property and humid air property database
TabularBackends.h
Go to the documentation of this file.
1 #ifndef TABULAR_BACKENDS_H
2 #define TABULAR_BACKENDS_H
3 
4 #include "AbstractState.h"
5 #include "CPmsgpack.h"
6 #include <msgpack/fbuffer.hpp>
8 #include "Exceptions.h"
9 #include "CoolProp.h"
10 #include <sstream>
11 #include "Configuration.h"
13 
18 #define LIST_OF_MATRICES \
19  X(T) \
20  X(p) \
21  X(rhomolar) \
22  X(hmolar) \
23  X(smolar) \
24  X(umolar) \
25  X(dTdx) \
26  X(dTdy) \
27  X(dpdx) \
28  X(dpdy) \
29  X(drhomolardx) \
30  X(drhomolardy) \
31  X(dhmolardx) \
32  X(dhmolardy) \
33  X(dsmolardx) \
34  X(dsmolardy) \
35  X(dumolardx) \
36  X(dumolardy) \
37  X(d2Tdx2) \
38  X(d2Tdxdy) \
39  X(d2Tdy2) \
40  X(d2pdx2) \
41  X(d2pdxdy) \
42  X(d2pdy2) \
43  X(d2rhomolardx2) \
44  X(d2rhomolardxdy) \
45  X(d2rhomolardy2) \
46  X(d2hmolardx2) \
47  X(d2hmolardxdy) \
48  X(d2hmolardy2) \
49  X(d2smolardx2) \
50  X(d2smolardxdy) \
51  X(d2smolardy2) \
52  X(d2umolardx2) \
53  X(d2umolardxdy) \
54  X(d2umolardy2) \
55  X(visc) \
56  X(cond)
57 
62 #define LIST_OF_SATURATION_VECTORS \
63  X(TL) \
64  X(pL) \
65  X(logpL) \
66  X(hmolarL) \
67  X(smolarL) \
68  X(umolarL) \
69  X(rhomolarL) \
70  X(logrhomolarL) \
71  X(viscL) \
72  X(condL) \
73  X(logviscL) \
74  X(TV) \
75  X(pV) \
76  X(logpV) \
77  X(hmolarV) \
78  X(smolarV) \
79  X(umolarV) \
80  X(rhomolarV) \
81  X(logrhomolarV) \
82  X(viscV) \
83  X(condV) \
84  X(logviscV) \
85  X(cpmolarV) \
86  X(cpmolarL) \
87  X(cvmolarV) \
88  X(cvmolarL) \
89  X(speed_soundL) \
90  X(speed_soundV)
91 
92 namespace CoolProp {
93 
95 {
96 
97  public:
98  int revision;
99 
101 
103 /* Use X macros to auto-generate the copying */
104 #define X(name) name = PED.name;
106 #undef X
107 /* Use X macros to auto-generate the copying */
108 #define X(name) name = PED.name;
110 #undef X
111  };
112 
113  std::map<std::string, std::vector<double>> vectors;
114  std::map<std::string, std::vector<std::vector<double>>> matrices;
115 
116  MSGPACK_DEFINE(revision, vectors, matrices); // write the member variables that you want to pack using msgpack
117 
119  void pack() {
120 /* Use X macros to auto-generate the packing code; each will look something like: matrices.insert(std::pair<std::string, std::vector<double> >("T", T)); */
121 #define X(name) vectors.insert(std::pair<std::string, std::vector<double>>(#name, name));
123 #undef X
124 /* Use X macros to auto-generate the packing code; each will look something like: matrices.insert(std::pair<std::string, std::vector<std::vector<CoolPropDbl> > >("T", T)); */
125 #define X(name) matrices.insert(std::pair<std::string, std::vector<std::vector<double>>>(#name, name));
127 #undef X
128  };
129  std::map<std::string, std::vector<double>>::iterator get_vector_iterator(const std::string& name) {
130  std::map<std::string, std::vector<double>>::iterator it = vectors.find(name);
131  if (it == vectors.end()) {
132  throw UnableToLoadError(format("could not find vector %s", name.c_str()));
133  }
134  return it;
135  }
136  std::map<std::string, std::vector<std::vector<double>>>::iterator get_matrix_iterator(const std::string& name) {
137  std::map<std::string, std::vector<std::vector<double>>>::iterator it = matrices.find(name);
138  if (it == matrices.end()) {
139  throw UnableToLoadError(format("could not find matrix %s", name.c_str()));
140  }
141  return it;
142  }
144  void unpack() {
145 /* Use X macros to auto-generate the unpacking code;
146  * each will look something like: T = get_vector_iterator("T")->second
147  */
148 #define X(name) name = get_vector_iterator(#name)->second;
150 #undef X
151 /* Use X macros to auto-generate the unpacking code;
152  * each will look something like: T = get_matrix_iterator("T")->second
153  **/
154 #define X(name) name = get_matrix_iterator(#name)->second;
156 #undef X
157  // Find the index of the point with the highest temperature
158  iTsat_max = std::distance(T.begin(), std::max_element(T.begin(), T.end()));
159  // Find the index of the point with the highest pressure
160  ipsat_max = std::distance(p.begin(), std::max_element(p.begin(), p.end()));
161  };
162  void deserialize(msgpack::object& deserialized) {
164  deserialized.convert(temp);
165  temp.unpack();
166  if (revision > temp.revision) {
167  throw ValueError(format("loaded revision [%d] is older than current revision [%d]", temp.revision, revision));
168  }
169  std::swap(*this, temp); // Swap if successful
170  };
171 };
172 
174 inline void mass_to_molar(parameters& param, double& conversion_factor, double molar_mass) {
175  conversion_factor = 1.0;
176  switch (param) {
177  case iDmass:
178  conversion_factor = molar_mass;
179  param = iDmolar;
180  break;
181  case iHmass:
182  conversion_factor /= molar_mass;
183  param = iHmolar;
184  break;
185  case iSmass:
186  conversion_factor /= molar_mass;
187  param = iSmolar;
188  break;
189  case iUmass:
190  conversion_factor /= molar_mass;
191  param = iUmolar;
192  break;
193  case iCvmass:
194  conversion_factor /= molar_mass;
195  param = iCvmolar;
196  break;
197  case iCpmass:
198  conversion_factor /= molar_mass;
199  param = iCpmolar;
200  break;
201  case iDmolar:
202  case iHmolar:
203  case iSmolar:
204  case iUmolar:
205  case iCvmolar:
206  case iCpmolar:
207  case iT:
208  case iP:
209  case ispeed_sound:
212  case iviscosity:
213  case iconductivity:
214  return;
215  default:
216  throw ValueError("TabularBackends::mass_to_molar - I don't know how to convert this parameter");
217  }
218 }
219 
225 {
226  public:
227  std::size_t N;
228  shared_ptr<CoolProp::AbstractState> AS;
229 
231  N = 1000;
232  revision = 1;
233  }
234 
236  void build(shared_ptr<CoolProp::AbstractState>& AS);
237 
238 /* Use X macros to auto-generate the variables; each will look something like: std::vector<double> T; */
239 #define X(name) std::vector<double> name;
241 #undef X
242 
243  int revision;
244  std::map<std::string, std::vector<double>> vectors;
245 
246  MSGPACK_DEFINE(revision, vectors); // write the member variables that you want to pack
247 
248  /***
249  * \brief Determine if a set of inputs are single-phase or inside the saturation table
250  * @param main The main variable that is being provided (currently T or P)
251  * @param mainval The value of the main variable that is being provided
252  * @param other The secondary variable
253  * @param val The value of the secondary variable
254  * @param iL The index associated with the nearest point for the liquid
255  * @param iV The index associated with the nearest point for the vapor
256  * @param yL The value associated with the nearest point for the liquid (based on interpolation)
257  * @param yV The value associated with the nearest point for the vapor (based on interpolation)
258 
259  \note If PQ or QT are inputs, yL and yV will correspond to the other main variable: p->T or T->p
260  */
261  bool is_inside(parameters main, double mainval, parameters other, double val, std::size_t& iL, std::size_t& iV, CoolPropDbl& yL,
262  CoolPropDbl& yV) {
263  std::vector<double>*yvecL = NULL, *yvecV = NULL;
264  switch (other) {
265  case iT:
266  yvecL = &TL;
267  yvecV = &TV;
268  break;
269  case iHmolar:
270  yvecL = &hmolarL;
271  yvecV = &hmolarV;
272  break;
273  case iQ:
274  yvecL = &TL;
275  yvecV = &TV;
276  break;
277  case iSmolar:
278  yvecL = &smolarL;
279  yvecV = &smolarV;
280  break;
281  case iUmolar:
282  yvecL = &umolarL;
283  yvecV = &umolarV;
284  break;
285  case iDmolar:
286  yvecL = &rhomolarL;
287  yvecV = &rhomolarV;
288  break;
289  default:
290  throw ValueError("invalid input for other in is_inside");
291  }
292 
293  // Trivial checks
294  if (main == iP) {
295  // If p is outside the range (ptriple, pcrit), considered to not be inside
296  double pmax = this->pV[pV.size() - 1], pmin = this->pV[0];
297  if (mainval > pmax || mainval < pmin) {
298  return false;
299  }
300  } else if (main == iT) {
301  // If T is outside the range (Tmin, Tcrit), considered to not be inside
302  double Tmax = this->TV[TV.size() - 1], Tmin = this->TV[0];
303  if (mainval > Tmax || mainval < Tmin) {
304  return false;
305  }
306  } else {
307  throw ValueError("invalid input for other in is_inside");
308  }
309 
310  // Now check based on a rough analysis using bounding pressure
311  std::size_t iLplus, iVplus;
312  // Find the indices (iL,iL+1) & (iV,iV+1) that bound the given pressure
313  // In general iV and iL will be the same, but if pseudo-pure, they might
314  // be different
315  if (main == iP) {
316  bisect_vector(pV, mainval, iV);
317  bisect_vector(pL, mainval, iL);
318  } else if (main == iT) {
319  bisect_vector(TV, mainval, iV);
320  bisect_vector(TL, mainval, iL);
321  } else {
322  throw ValueError(format("For now, main input in is_inside must be T or p"));
323  }
324 
325  iVplus = std::min(iV + 1, N - 1);
326  iLplus = std::min(iL + 1, N - 1);
327  if (other == iQ) {
328  // Actually do "saturation" call using cubic interpolation
329  if (iVplus < 3) {
330  iVplus = 3;
331  }
332  if (iLplus < 3) {
333  iLplus = 3;
334  }
335  if (main == iP) {
336  double logp = log(mainval);
337  // Calculate temperature
338  yV = CubicInterp(logpV, TV, iVplus - 3, iVplus - 2, iVplus - 1, iVplus, logp);
339  yL = CubicInterp(logpL, TL, iLplus - 3, iLplus - 2, iLplus - 1, iLplus, logp);
340  } else if (main == iT) {
341  // Calculate pressure
342  yV = exp(CubicInterp(TV, logpV, iVplus - 3, iVplus - 2, iVplus - 1, iVplus, mainval));
343  yL = exp(CubicInterp(TL, logpL, iLplus - 3, iLplus - 2, iLplus - 1, iLplus, mainval));
344  }
345  return true;
346  }
347  // Find the bounding values for the other variable
348  double ymin = min4((*yvecL)[iL], (*yvecL)[iLplus], (*yvecV)[iV], (*yvecV)[iVplus]);
349  double ymax = max4((*yvecL)[iL], (*yvecL)[iLplus], (*yvecV)[iV], (*yvecV)[iVplus]);
350  if (val < ymin || val > ymax) {
351  return false;
352  }
353  // Actually do "saturation" call using cubic interpolation
354  if (iVplus < 3) {
355  iVplus = 3;
356  }
357  if (iLplus < 3) {
358  iLplus = 3;
359  }
360  if (main == iP) {
361  double logp = log(mainval);
362  yV = CubicInterp(logpV, *yvecV, iVplus - 3, iVplus - 2, iVplus - 1, iVplus, logp);
363  yL = CubicInterp(logpL, *yvecL, iLplus - 3, iLplus - 2, iLplus - 1, iLplus, logp);
364  } else if (main == iT) {
365  yV = CubicInterp(TV, *yvecV, iVplus - 3, iVplus - 2, iVplus - 1, iVplus, mainval);
366  yL = CubicInterp(TL, *yvecL, iLplus - 3, iLplus - 2, iLplus - 1, iLplus, mainval);
367  }
368 
369  if (!is_in_closed_range(yV, yL, static_cast<CoolPropDbl>(val))) {
370  return false;
371  } else {
372  iL = iLplus - 1;
373  iV = iVplus - 1;
374  return true;
375  }
376  }
378  void resize(std::size_t N) {
379 /* Use X macros to auto-generate the code; each will look something like: T.resize(N); std::fill(T.begin(), T.end(), _HUGE); */
380 #define X(name) \
381  name.resize(N); \
382  std::fill(name.begin(), name.end(), _HUGE);
384 #undef X
385  };
387  void pack() {
388 /* Use X macros to auto-generate the packing code; each will look something like: matrices.insert(std::pair<std::vector<std::vector<double> > >("T", T)); */
389 #define X(name) vectors.insert(std::pair<std::string, std::vector<double>>(#name, name));
391 #undef X
392  };
393  std::map<std::string, std::vector<double>>::iterator get_vector_iterator(const std::string& name) {
394  std::map<std::string, std::vector<double>>::iterator it = vectors.find(name);
395  if (it == vectors.end()) {
396  throw UnableToLoadError(format("could not find vector %s", name.c_str()));
397  }
398  return it;
399  }
401  void unpack() {
402 /* Use X macros to auto-generate the unpacking code; each will look something like: T = get_vector_iterator("T")->second */
403 #define X(name) name = get_vector_iterator(#name)->second;
405 #undef X
406  N = TL.size();
407  };
408  void deserialize(msgpack::object& deserialized) {
410  deserialized.convert(temp);
411  temp.unpack();
412  if (N != temp.N) {
413  throw ValueError(format("old [%d] and new [%d] sizes don't agree", temp.N, N));
414  } else if (revision > temp.revision) {
415  throw ValueError(format("loaded revision [%d] is older than current revision [%d]", temp.revision, revision));
416  }
417  std::swap(*this, temp); // Swap
418  this->AS = temp.AS; // Reconnect the AbstractState pointer
419  };
420  double evaluate(parameters output, double p_or_T, double Q, std::size_t iL, std::size_t iV) {
421  if (iL <= 2) {
422  iL = 2;
423  } else if (iL + 1 == N) {
424  iL = N - 2;
425  }
426  if (iV <= 2) {
427  iV = 2;
428  } else if (iV + 1 == N) {
429  iV = N - 2;
430  }
431  double logp = log(p_or_T);
432  switch (output) {
433  case iP: {
434  double _logpV = CubicInterp(this->TV, logpV, iV - 2, iV - 1, iV, iV + 1, p_or_T);
435  double _logpL = CubicInterp(this->TL, logpL, iL - 2, iL - 1, iL, iL + 1, p_or_T);
436  return Q * exp(_logpV) + (1 - Q) * exp(_logpL);
437  }
438  case iT: {
439  double TV = CubicInterp(logpV, this->TV, iV - 2, iV - 1, iV, iV + 1, logp);
440  double TL = CubicInterp(logpL, this->TL, iL - 2, iL - 1, iL, iL + 1, logp);
441  return Q * TV + (1 - Q) * TL;
442  }
443  case iSmolar: {
444  double sV = CubicInterp(logpV, smolarV, iV - 2, iV - 1, iV, iV + 1, logp);
445  double sL = CubicInterp(logpL, smolarL, iL - 2, iL - 1, iL, iL + 1, logp);
446  return Q * sV + (1 - Q) * sL;
447  }
448  case iHmolar: {
449  double hV = CubicInterp(logpV, hmolarV, iV - 2, iV - 1, iV, iV + 1, logp);
450  double hL = CubicInterp(logpL, hmolarL, iL - 2, iL - 1, iL, iL + 1, logp);
451  return Q * hV + (1 - Q) * hL;
452  }
453  case iUmolar: {
454  double uV = CubicInterp(logpV, umolarV, iV - 2, iV - 1, iV, iV + 1, logp);
455  double uL = CubicInterp(logpL, umolarL, iL - 2, iL - 1, iL, iL + 1, logp);
456  return Q * uV + (1 - Q) * uL;
457  }
458  case iDmolar: {
459  double rhoV = exp(CubicInterp(logpV, logrhomolarV, iV - 2, iV - 1, iV, iV + 1, logp));
460  double rhoL = exp(CubicInterp(logpL, logrhomolarL, iL - 2, iL - 1, iL, iL + 1, logp));
461  if (!ValidNumber(rhoV)) {
462  throw ValueError("rhoV is invalid");
463  }
464  if (!ValidNumber(rhoL)) {
465  throw ValueError("rhoL is invalid");
466  }
467  return 1 / (Q / rhoV + (1 - Q) / rhoL);
468  }
469  case iconductivity: {
470  double kV = CubicInterp(logpV, condV, iV - 2, iV - 1, iV, iV + 1, logp);
471  double kL = CubicInterp(logpL, condL, iL - 2, iL - 1, iL, iL + 1, logp);
472  if (!ValidNumber(kV)) {
473  throw ValueError("kV is invalid");
474  }
475  if (!ValidNumber(kL)) {
476  throw ValueError("kL is invalid");
477  }
478  return Q * kV + (1 - Q) * kL;
479  }
480  case iviscosity: {
481  double muV = exp(CubicInterp(logpV, logviscV, iV - 2, iV - 1, iV, iV + 1, logp));
482  double muL = exp(CubicInterp(logpL, logviscL, iL - 2, iL - 1, iL, iL + 1, logp));
483  if (!ValidNumber(muV)) {
484  throw ValueError("muV is invalid");
485  }
486  if (!ValidNumber(muL)) {
487  throw ValueError("muL is invalid");
488  }
489  return 1 / (Q / muV + (1 - Q) / muL);
490  }
491  case iCpmolar: {
492  double cpV = CubicInterp(logpV, cpmolarV, iV - 2, iV - 1, iV, iV + 1, logp);
493  double cpL = CubicInterp(logpL, cpmolarL, iL - 2, iL - 1, iL, iL + 1, logp);
494  if (!ValidNumber(cpV)) {
495  throw ValueError("cpV is invalid");
496  }
497  if (!ValidNumber(cpL)) {
498  throw ValueError("cpL is invalid");
499  }
500  return Q * cpV + (1 - Q) * cpL;
501  }
502  case iCvmolar: {
503  double cvV = CubicInterp(logpV, cvmolarV, iV - 2, iV - 1, iV, iV + 1, logp);
504  double cvL = CubicInterp(logpL, cvmolarL, iL - 2, iL - 1, iL, iL + 1, logp);
505  if (!ValidNumber(cvV)) {
506  throw ValueError("cvV is invalid");
507  }
508  if (!ValidNumber(cvL)) {
509  throw ValueError("cvL is invalid");
510  }
511  return Q * cvV + (1 - Q) * cvL;
512  }
513  case ispeed_sound: {
514  double wV = CubicInterp(logpV, speed_soundV, iV - 2, iV - 1, iV, iV + 1, logp);
515  double wL = CubicInterp(logpL, speed_soundL, iL - 2, iL - 1, iL, iL + 1, logp);
516  if (!ValidNumber(wV)) {
517  throw ValueError("wV is invalid");
518  }
519  if (!ValidNumber(wL)) {
520  throw ValueError("wL is invalid");
521  }
522  return Q * wV + (1 - Q) * wL;
523  }
524  default:
525  throw ValueError("Output variable for evaluate is invalid");
526  }
527  };
536  double first_saturation_deriv(parameters Of1, parameters Wrt1, int Q, double val, std::size_t i) {
537  if (i < 2 || i > TL.size() - 2) {
538  throw ValueError(format("Invalid index (%d) to calc_first_saturation_deriv in TabularBackends", i));
539  }
540  std::vector<double>*x, *y;
541  // Connect pointers for each vector
542  switch (Wrt1) {
543  case iT:
544  x = (Q == 0) ? &TL : &TV;
545  break;
546  case iP:
547  x = (Q == 0) ? &pL : &pV;
548  break;
549  default:
550  throw ValueError(format("Key for Wrt1 is invalid in calc_first_saturation_deriv"));
551  }
552  CoolPropDbl factor = 1.0;
553  switch (Of1) {
554  case iT:
555  y = (Q == 0) ? &TL : &TV;
556  break;
557  case iP:
558  y = (Q == 0) ? &pL : &pV;
559  break;
560  case iDmolar:
561  y = (Q == 0) ? &rhomolarL : &rhomolarV;
562  break;
563  case iHmolar:
564  y = (Q == 0) ? &hmolarL : &hmolarV;
565  break;
566  case iSmolar:
567  y = (Q == 0) ? &smolarL : &smolarV;
568  break;
569  case iUmolar:
570  y = (Q == 0) ? &umolarL : &umolarV;
571  break;
572  case iDmass:
573  y = (Q == 0) ? &rhomolarL : &rhomolarV;
574  factor = AS->molar_mass();
575  break;
576  case iHmass:
577  y = (Q == 0) ? &hmolarL : &hmolarV;
578  factor = 1 / AS->molar_mass();
579  break;
580  case iSmass:
581  y = (Q == 0) ? &smolarL : &smolarV;
582  factor = 1 / AS->molar_mass();
583  break;
584  case iUmass:
585  y = (Q == 0) ? &umolarL : &umolarV;
586  factor = 1 / AS->molar_mass();
587  break;
588  default:
589  throw ValueError(format("Key for Of1 is invalid in calc_first_saturation_deriv"));
590  }
591  return CubicInterpFirstDeriv((*x)[i - 2], (*x)[i - 1], (*x)[i], (*x)[i + 1], (*y)[i - 2], (*y)[i - 1], (*y)[i], (*y)[i + 1], val) * factor;
592  };
593  //calc_first_two_phase_deriv(parameters Of, parameters Wrt, parameters Constant);
594 };
595 
601 {
602 
603  public:
604  std::size_t Nx, Ny;
606  shared_ptr<CoolProp::AbstractState> AS;
607  std::vector<double> xvec, yvec;
608  std::vector<std::vector<std::size_t>> nearest_neighbor_i, nearest_neighbor_j;
609  bool logx, logy;
610  double xmin, ymin, xmax, ymax;
611 
612  virtual void set_limits() = 0;
613 
615  Nx = 200;
616  Ny = 200;
617  revision = 0;
620  logx = false;
621  logy = false;
622  xmin = _HUGE;
623  xmax = _HUGE;
624  ymin = _HUGE;
625  ymax = _HUGE;
626  }
627 
628 /* Use X macros to auto-generate the variables; each will look something like: std::vector< std::vector<double> > T; */
629 #define X(name) std::vector<std::vector<double>> name;
631 #undef X
632  int revision;
633  std::map<std::string, std::vector<std::vector<double>>> matrices;
635  void build(shared_ptr<CoolProp::AbstractState>& AS);
636 
637  MSGPACK_DEFINE(revision, matrices, xmin, xmax, ymin, ymax); // write the member variables that you want to pack
639  void resize(std::size_t Nx, std::size_t Ny) {
640 /* Use X macros to auto-generate the code; each will look something like: T.resize(Nx, std::vector<double>(Ny, _HUGE)); */
641 #define X(name) name.resize(Nx, std::vector<double>(Ny, _HUGE));
643 #undef X
645  };
647  void make_axis_vectors(void) {
648  if (logx) {
649  xvec = logspace(xmin, xmax, Nx);
650  } else {
651  xvec = linspace(xmin, xmax, Nx);
652  }
653  if (logy) {
654  yvec = logspace(ymin, ymax, Ny);
655  } else {
656  yvec = linspace(ymin, ymax, Ny);
657  }
658  };
660  void make_good_neighbors(void) {
661  nearest_neighbor_i.resize(Nx, std::vector<std::size_t>(Ny, std::numeric_limits<std::size_t>::max()));
662  nearest_neighbor_j.resize(Nx, std::vector<std::size_t>(Ny, std::numeric_limits<std::size_t>::max()));
663  for (std::size_t i = 0; i < xvec.size(); ++i) {
664  for (std::size_t j = 0; j < yvec.size(); ++j) {
665  nearest_neighbor_i[i][j] = i;
666  nearest_neighbor_j[i][j] = j;
667  if (!ValidNumber(T[i][j])) {
668  int xoffsets[] = {-1, 1, 0, 0, -1, 1, 1, -1};
669  int yoffsets[] = {0, 0, 1, -1, -1, -1, 1, 1};
670  // Length of offset
671  std::size_t N = sizeof(xoffsets) / sizeof(xoffsets[0]);
672  for (std::size_t k = 0; k < N; ++k) {
673  std::size_t iplus = i + xoffsets[k];
674  std::size_t jplus = j + yoffsets[k];
675  if (0 < iplus && iplus < Nx - 1 && 0 < jplus && jplus < Ny - 1 && ValidNumber(T[iplus][jplus])) {
676  nearest_neighbor_i[i][j] = iplus;
677  nearest_neighbor_j[i][j] = jplus;
678  break;
679  }
680  }
681  }
682  }
683  }
684  };
686  void pack() {
687 /* Use X macros to auto-generate the packing code; each will look something like: matrices.insert(std::pair<std::vector<std::vector<double> > >("T", T)); */
688 #define X(name) matrices.insert(std::pair<std::string, std::vector<std::vector<double>>>(#name, name));
690 #undef X
691  };
692  std::map<std::string, std::vector<std::vector<double>>>::iterator get_matrices_iterator(const std::string& name) {
693  std::map<std::string, std::vector<std::vector<double>>>::iterator it = matrices.find(name);
694  if (it == matrices.end()) {
695  throw UnableToLoadError(format("could not find matrix %s", name.c_str()));
696  }
697  return it;
698  }
700  void unpack() {
701 /* Use X macros to auto-generate the unpacking code; each will look something like: T = *(matrices.find("T")).second */
702 #define X(name) name = get_matrices_iterator(#name)->second;
704 #undef X
705  Nx = T.size();
706  Ny = T[0].size();
709  };
711  bool native_inputs_are_in_range(double x, double y) {
712  double e = 10 * DBL_EPSILON;
713  return x >= xmin - e && x <= xmax + e && y >= ymin - e && y <= ymax + e;
714  }
718  void find_native_nearest_neighbor(double x, double y, std::size_t& i, std::size_t& j) {
719  bisect_vector(xvec, x, i);
720  if (i != Nx - 1) {
721  if (!logx) {
722  if (x > (xvec[i] + xvec[i + 1]) / 2.0) {
723  i++;
724  }
725  } else {
726  if (x > sqrt(xvec[i] * xvec[i + 1])) {
727  i++;
728  }
729  }
730  }
731  bisect_vector(yvec, y, j);
732  if (j != Ny - 1) {
733  if (!logy) {
734  if (y > (yvec[j] + yvec[j + 1]) / 2.0) {
735  j++;
736  }
737  } else {
738  if (y > sqrt(yvec[j] * yvec[j + 1])) {
739  j++;
740  }
741  }
742  }
743  }
745  void find_nearest_neighbor(parameters givenkey, double givenval, parameters otherkey, double otherval, std::size_t& i, std::size_t& j) {
746  if (givenkey == ykey) {
747  bisect_vector(yvec, givenval, j);
748  // This one is problematic because we need to make a slice against the grain in the "matrix"
749  // which requires a slightly different algorithm
750  try {
751  bisect_segmented_vector_slice(get(otherkey), j, otherval, i);
752  } catch (...) {
753  // Now we go for a less intelligent solution, we simply try to find the one that is the closest
754  const std::vector<std::vector<double>>& mat = get(otherkey);
755  double closest_diff = 1e20;
756  std::size_t closest_i = 0;
757  for (std::size_t index = 0; index < mat.size(); ++index) {
758  double diff = std::abs(mat[index][j] - otherval);
759  if (diff < closest_diff) {
760  closest_diff = diff;
761  closest_i = index;
762  }
763  }
764  i = closest_i;
765  }
766  } else if (givenkey == xkey) {
767  bisect_vector(xvec, givenval, i);
768  // This one is fine because we now end up with a vector<double> in the other variable
769  const std::vector<std::vector<double>>& v = get(otherkey);
770  bisect_vector(v[i], otherval, j);
771  }
772  }
775  void find_native_nearest_good_neighbor(double x, double y, std::size_t& i, std::size_t& j) {
776  // Get the node that is closest
777  find_native_nearest_neighbor(x, y, i, j);
778  // Check whether found node is good
779  if (!ValidNumber(T[i][j])) {
780  // If not, find its nearest good neighbor
781  // (nearest good neighbors are precalculated and cached)
782  std::size_t inew = nearest_neighbor_i[i][j];
783  std::size_t jnew = nearest_neighbor_j[i][j];
784  i = inew;
785  j = jnew;
786  }
787  }
790  void find_native_nearest_good_cell(double x, double y, std::size_t& i, std::size_t& j) {
791  bisect_vector(xvec, x, i);
792  bisect_vector(yvec, y, j);
793  }
794  const std::vector<std::vector<double>>& get(parameters key) {
795  switch (key) {
796  case iDmolar:
797  return rhomolar;
798  case iT:
799  return T;
800  case iUmolar:
801  return umolar;
802  case iHmolar:
803  return hmolar;
804  case iSmolar:
805  return smolar;
806  case iP:
807  return p;
808  case iviscosity:
809  return visc;
810  case iconductivity:
811  return cond;
812  default:
813  throw KeyError(format("invalid key"));
814  }
815  }
816 };
817 
820 {
821  public:
823  xkey = iHmolar;
824  ykey = iP;
825  logy = true;
826  logx = false;
827  };
828  void set_limits() {
829  if (this->AS.get() == NULL) {
830  throw ValueError("AS is not yet set");
831  }
832  CoolPropDbl Tmin = std::max(AS->Ttriple(), AS->Tmin());
833  // Minimum enthalpy is the saturated liquid enthalpy
834  AS->update(QT_INPUTS, 0, Tmin);
835  xmin = AS->hmolar();
836  ymin = AS->p();
837 
838  // Check both the enthalpies at the Tmax isotherm to see whether to use low or high pressure
839  AS->update(DmolarT_INPUTS, 1e-10, 1.499 * AS->Tmax());
840  CoolPropDbl xmax1 = AS->hmolar();
841  AS->update(PT_INPUTS, AS->pmax(), 1.499 * AS->Tmax());
842  CoolPropDbl xmax2 = AS->hmolar();
843  xmax = std::max(xmax1, xmax2);
844 
845  ymax = AS->pmax();
846  }
847  void deserialize(msgpack::object& deserialized) {
848  LogPHTable temp;
849  deserialized.convert(temp);
850  temp.unpack();
851  if (Nx != temp.Nx || Ny != temp.Ny) {
852  throw ValueError(format("old [%dx%d] and new [%dx%d] dimensions don't agree", temp.Nx, temp.Ny, Nx, Ny));
853  } else if (revision > temp.revision) {
854  throw ValueError(format("loaded revision [%d] is older than current revision [%d]", temp.revision, revision));
855  } else if ((std::abs(xmin) > 1e-10 && std::abs(xmax) > 1e-10)
856  && (std::abs(temp.xmin - xmin) / xmin > 1e-6 || std::abs(temp.xmax - xmax) / xmax > 1e-6)) {
857  throw ValueError(format("Current limits for x [%g,%g] do not agree with loaded limits [%g,%g]", xmin, xmax, temp.xmin, temp.xmax));
858  } else if ((std::abs(ymin) > 1e-10 && std::abs(ymax) > 1e-10)
859  && (std::abs(temp.ymin - ymin) / ymin > 1e-6 || std::abs(temp.ymax - ymax) / ymax > 1e-6)) {
860  throw ValueError(format("Current limits for y [%g,%g] do not agree with loaded limits [%g,%g]", ymin, ymax, temp.ymin, temp.ymax));
861  }
862  std::swap(*this, temp); // Swap
863  this->AS = temp.AS; // Reconnect the AbstractState pointer
864  };
865 };
868 {
869  public:
871  xkey = iT;
872  ykey = iP;
873  logy = true;
874  logx = false;
875  xmin = _HUGE;
876  ymin = _HUGE;
877  xmax = _HUGE;
878  ymax = _HUGE;
879  };
880  void set_limits() {
881  if (this->AS.get() == NULL) {
882  throw ValueError("AS is not yet set");
883  }
884  CoolPropDbl Tmin = std::max(AS->Ttriple(), AS->Tmin());
885  AS->update(QT_INPUTS, 0, Tmin);
886  xmin = Tmin;
887  ymin = AS->p();
888 
889  xmax = AS->Tmax() * 1.499;
890  ymax = AS->pmax();
891  }
892  void deserialize(msgpack::object& deserialized) {
893  LogPTTable temp;
894  deserialized.convert(temp);
895  temp.unpack();
896  if (Nx != temp.Nx || Ny != temp.Ny) {
897  throw ValueError(format("old [%dx%d] and new [%dx%d] dimensions don't agree", temp.Nx, temp.Ny, Nx, Ny));
898  } else if (revision > temp.revision) {
899  throw ValueError(format("loaded revision [%d] is older than current revision [%d]", temp.revision, revision));
900  } else if ((std::abs(xmin) > 1e-10 && std::abs(xmax) > 1e-10)
901  && (std::abs(temp.xmin - xmin) / xmin > 1e-6 || std::abs(temp.xmax - xmax) / xmax > 1e-6)) {
902  throw ValueError(format("Current limits for x [%g,%g] do not agree with loaded limits [%g,%g]", xmin, xmax, temp.xmin, temp.xmax));
903  } else if ((std::abs(ymin) > 1e-10 && std::abs(ymax) > 1e-10)
904  && (std::abs(temp.ymin - ymin) / ymin > 1e-6 || std::abs(temp.ymax - ymax) / ymax > 1e-6)) {
905  throw ValueError(format("Current limits for y [%g,%g] do not agree with loaded limits [%g,%g]", ymin, ymax, temp.ymin, temp.ymax));
906  }
907  std::swap(*this, temp); // Swap
908  this->AS = temp.AS; // Reconnect the AbstractState pointer
909  };
910 };
911 
915 {
916  private:
917  std::size_t alt_i, alt_j;
918  bool _valid, _has_valid_neighbor;
919 
920  public:
923  _valid = false;
924  _has_valid_neighbor = false;
925  dx_dxhat = _HUGE;
926  dy_dyhat = _HUGE;
927  alt_i = 9999999;
928  alt_j = 9999999;
929  }
930  std::vector<double> T, rhomolar, hmolar, p, smolar, umolar;
932  const std::vector<double>& get(const parameters params) const {
933  switch (params) {
934  case iT:
935  return T;
936  case iP:
937  return p;
938  case iDmolar:
939  return rhomolar;
940  case iHmolar:
941  return hmolar;
942  case iSmolar:
943  return smolar;
944  case iUmolar:
945  return umolar;
946  default:
947  throw KeyError(format("Invalid key to get() function of CellCoeffs"));
948  }
949  };
951  void set(parameters params, const std::vector<double>& mat) {
952  switch (params) {
953  case iT:
954  T = mat;
955  break;
956  case iP:
957  p = mat;
958  break;
959  case iDmolar:
960  rhomolar = mat;
961  break;
962  case iHmolar:
963  hmolar = mat;
964  break;
965  case iSmolar:
966  smolar = mat;
967  break;
968  case iUmolar:
969  umolar = mat;
970  break;
971  default:
972  throw KeyError(format("Invalid key to set() function of CellCoeffs"));
973  }
974  };
976  bool valid() const {
977  return _valid;
978  };
980  void set_valid() {
981  _valid = true;
982  };
984  void set_invalid() {
985  _valid = false;
986  };
988  void set_alternate(std::size_t i, std::size_t j) {
989  alt_i = i;
990  alt_j = j;
991  _has_valid_neighbor = true;
992  }
994  void get_alternate(std::size_t& i, std::size_t& j) const {
995  if (_has_valid_neighbor) {
996  i = alt_i;
997  j = alt_j;
998  } else {
999  throw ValueError("No valid neighbor");
1000  }
1001  }
1003  bool has_valid_neighbor() const {
1004  return _has_valid_neighbor;
1005  }
1006 };
1007 
1010 {
1011  public:
1017  std::vector<std::vector<CellCoeffs>> coeffs_ph, coeffs_pT;
1018 
1020  tables_loaded = false;
1021  }
1023  void write_tables(const std::string& path_to_tables);
1025  void load_tables(const std::string& path_to_tables, shared_ptr<CoolProp::AbstractState>& AS);
1027  void build_tables(shared_ptr<CoolProp::AbstractState>& AS);
1029  void build_coeffs(SinglePhaseGriddedTableData& table, std::vector<std::vector<CellCoeffs>>& coeffs);
1030 };
1031 
1033 {
1034  private:
1035  std::map<std::string, TabularDataSet> data;
1036 
1037  public:
1039  std::string path_to_tables(shared_ptr<CoolProp::AbstractState>& AS) {
1040  std::vector<std::string> fluids = AS->fluid_names();
1041  std::vector<CoolPropDbl> fractions = AS->get_mole_fractions();
1042  std::vector<std::string> components;
1043  for (std::size_t i = 0; i < fluids.size(); ++i) {
1044  components.push_back(format("%s[%0.10Lf]", fluids[i].c_str(), fractions[i]));
1045  }
1046  std::string table_directory = get_home_dir() + "/.CoolProp/Tables/";
1047  std::string alt_table_directory = get_config_string(ALTERNATIVE_TABLES_DIRECTORY);
1048  if (!alt_table_directory.empty()) {
1049  table_directory = alt_table_directory;
1050  }
1051  return table_directory + AS->backend_name() + "(" + strjoin(components, "&") + ")";
1052  }
1054  TabularDataSet* get_set_of_tables(shared_ptr<AbstractState>& AS, bool& loaded);
1055 };
1056 
1064 {
1065  protected:
1069  {
1073  };
1077  std::vector<std::vector<double>> const* z;
1078  std::vector<std::vector<double>> const* dzdx;
1079  std::vector<std::vector<double>> const* dzdy;
1080  std::vector<std::vector<double>> const* d2zdx2;
1081  std::vector<std::vector<double>> const* d2zdxdy;
1082  std::vector<std::vector<double>> const* d2zdy2;
1083  std::vector<CoolPropDbl> mole_fractions;
1084 
1085  public:
1086  shared_ptr<CoolProp::AbstractState> AS;
1087  TabularBackend(shared_ptr<CoolProp::AbstractState> AS) : tables_loaded(false), using_single_phase_table(false), is_mixture(false), AS(AS) {
1089  // Flush the cached indices (set to large number)
1090  cached_single_phase_i = std::numeric_limits<std::size_t>::max();
1091  cached_single_phase_j = std::numeric_limits<std::size_t>::max();
1092  cached_saturation_iL = std::numeric_limits<std::size_t>::max();
1093  cached_saturation_iV = std::numeric_limits<std::size_t>::max();
1094  z = NULL;
1095  dzdx = NULL;
1096  dzdy = NULL;
1097  d2zdx2 = NULL;
1098  d2zdxdy = NULL;
1099  d2zdy2 = NULL;
1100  dataset = NULL;
1102  };
1103 
1104  // None of the tabular methods are available from the high-level interface
1106  return false;
1107  }
1108 
1109  std::string calc_name(void) {
1110  return AS->name();
1111  }
1112  std::vector<std::string> calc_fluid_names(void) {
1113  return AS->fluid_names();
1114  }
1115 
1117  // Connect the pointers based on the output variable desired
1118  switch (output) {
1119  case iT:
1120  z = &table.T;
1121  dzdx = &table.dTdx;
1122  dzdy = &table.dTdy;
1123  d2zdxdy = &table.d2Tdxdy;
1124  d2zdx2 = &table.d2Tdx2;
1125  d2zdy2 = &table.d2Tdy2;
1126  break;
1127  case iDmolar:
1128  z = &table.rhomolar;
1129  dzdx = &table.drhomolardx;
1130  dzdy = &table.drhomolardy;
1131  d2zdxdy = &table.d2rhomolardxdy;
1132  d2zdx2 = &table.d2rhomolardx2;
1133  d2zdy2 = &table.d2rhomolardy2;
1134  break;
1135  case iSmolar:
1136  z = &table.smolar;
1137  dzdx = &table.dsmolardx;
1138  dzdy = &table.dsmolardy;
1139  d2zdxdy = &table.d2smolardxdy;
1140  d2zdx2 = &table.d2smolardx2;
1141  d2zdy2 = &table.d2smolardy2;
1142  break;
1143  case iHmolar:
1144  z = &table.hmolar;
1145  dzdx = &table.dhmolardx;
1146  dzdy = &table.dhmolardy;
1147  d2zdxdy = &table.d2hmolardxdy;
1148  d2zdx2 = &table.d2hmolardx2;
1149  d2zdy2 = &table.d2hmolardy2;
1150  break;
1151  case iUmolar:
1152  z = &table.umolar;
1153  dzdx = &table.dumolardx;
1154  dzdy = &table.dumolardy;
1155  d2zdxdy = &table.d2umolardxdy;
1156  d2zdx2 = &table.d2umolardx2;
1157  d2zdy2 = &table.d2umolardy2;
1158  break;
1159  case iviscosity:
1160  z = &table.visc;
1161  break;
1162  case iconductivity:
1163  z = &table.cond;
1164  break;
1165  default:
1166  throw ValueError();
1167  }
1168  }
1170 
1172  if (p() > p_critical()) {
1173  if (T() > T_critical()) {
1175  } else {
1177  }
1178  } else {
1179  if (T() > T_critical()) {
1181  } else {
1182  // Liquid or vapor
1183  if (rhomolar() > rhomolar_critical()) {
1185  } else {
1186  _phase = iphase_gas;
1187  }
1188  }
1189  }
1190  }
1195  void calc_specify_phase(phases phase_index) {
1196  imposed_phase_index = phase_index;
1197  };
1198 
1203  };
1204 
1205  virtual double evaluate_single_phase_phmolar(parameters output, std::size_t i, std::size_t j) = 0;
1206  virtual double evaluate_single_phase_pT(parameters output, std::size_t i, std::size_t j) = 0;
1207  virtual double evaluate_single_phase_phmolar_transport(parameters output, std::size_t i, std::size_t j) = 0;
1208  virtual double evaluate_single_phase_pT_transport(parameters output, std::size_t i, std::size_t j) = 0;
1209  virtual double evaluate_single_phase_phmolar_derivative(parameters output, std::size_t i, std::size_t j, std::size_t Nx, std::size_t Ny) = 0;
1210  virtual double evaluate_single_phase_pT_derivative(parameters output, std::size_t i, std::size_t j, std::size_t Nx, std::size_t Ny) = 0;
1211 
1213  virtual void find_native_nearest_good_indices(SinglePhaseGriddedTableData& table, const std::vector<std::vector<CellCoeffs>>& coeffs, double x,
1214  double y, std::size_t& i, std::size_t& j) = 0;
1216  virtual void find_nearest_neighbor(SinglePhaseGriddedTableData& table, const std::vector<std::vector<CellCoeffs>>& coeffs,
1217  const parameters variable1, const double value1, const parameters other, const double otherval, std::size_t& i,
1218  std::size_t& j) = 0;
1220  virtual void invert_single_phase_x(const SinglePhaseGriddedTableData& table, const std::vector<std::vector<CellCoeffs>>& coeffs,
1221  parameters output, double x, double y, std::size_t i, std::size_t j) = 0;
1223  virtual void invert_single_phase_y(const SinglePhaseGriddedTableData& table, const std::vector<std::vector<CellCoeffs>>& coeffs,
1224  parameters output, double x, double y, std::size_t i, std::size_t j) = 0;
1225 
1227  return _phase;
1228  }
1230  return this->AS->T_critical();
1231  };
1233  return this->AS->Ttriple();
1234  };
1236  return this->AS->p_triple();
1237  };
1239  return this->AS->pmax();
1240  };
1242  return this->AS->Tmax();
1243  };
1245  return this->AS->Tmin();
1246  };
1248  return this->AS->p_critical();
1249  }
1251  return this->AS->rhomolar_critical();
1252  }
1254  return true;
1255  }
1257  return false;
1258  }
1260  return false;
1261  }
1262  void update(CoolProp::input_pairs input_pair, double Value1, double Value2);
1263  void set_mole_fractions(const std::vector<CoolPropDbl>& mole_fractions) {
1264  this->AS->set_mole_fractions(mole_fractions);
1265  };
1266  void set_mass_fractions(const std::vector<CoolPropDbl>& mass_fractions) {
1267  throw NotImplementedError("set_mass_fractions not implemented for Tabular backends");
1268  };
1269  const std::vector<CoolPropDbl>& get_mole_fractions() {
1270  return AS->get_mole_fractions();
1271  };
1272  const std::vector<CoolPropDbl> calc_mass_fractions(void) {
1273  return AS->get_mass_fractions();
1274  };
1275 
1277  return AS->molar_mass();
1278  };
1279 
1282 
1284  std::string path_to_tables(void);
1286  void load_tables();
1287  void pack_matrices() {
1292  single_phase_logph.pack();
1293  single_phase_logpT.pack();
1294  pure_saturation.pack();
1295  phase_envelope.pack();
1296  }
1298  void write_tables();
1299 
1300  CoolPropDbl phase_envelope_sat(const PhaseEnvelopeData& env, parameters output, parameters iInput1, double value1) {
1301  const PhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
1302  CoolPropDbl yL = PhaseEnvelopeRoutines::evaluate(phase_envelope, output, iInput1, value1, cached_saturation_iL);
1303  CoolPropDbl yV = PhaseEnvelopeRoutines::evaluate(phase_envelope, output, iInput1, value1, cached_saturation_iV);
1304  return _Q * yV + (1 - _Q) * yL;
1305  }
1307  this->AS->set_T(_T);
1308  return this->AS->cp0molar();
1309  }
1312  this->AS->set_T(_T);
1313  return this->AS->surface_tension();
1314  this->AS->set_T(_HUGE);
1315  }
1316  CoolPropDbl calc_p(void);
1317  CoolPropDbl calc_T(void);
1318  CoolPropDbl calc_rhomolar(void);
1319  CoolPropDbl calc_hmolar(void);
1320  CoolPropDbl calc_smolar(void);
1321  CoolPropDbl calc_umolar(void);
1322  CoolPropDbl calc_cpmolar(void);
1323  CoolPropDbl calc_cvmolar(void);
1333 
1336 
1337  void check_tables() {
1338  if (!tables_loaded) {
1339  try {
1341  load_tables();
1342  // Set the flag saying tables have been successfully loaded
1343  tables_loaded = true;
1344  } catch (CoolProp::UnableToLoadError& e) {
1345  if (get_debug_level() > 0) {
1346  std::cout << format("Table loading failed with error: %s\n", e.what());
1347  }
1349  std::string table_path = path_to_tables();
1350 #if defined(__ISWINDOWS__)
1351  double directory_size_in_GB = CalculateDirSize(std::wstring(table_path.begin(), table_path.end())) / POW3(1024.0);
1352 #else
1353  double directory_size_in_GB = CalculateDirSize(table_path) / POW3(1024.0);
1354 #endif
1355  double allowed_size_in_GB = get_config_double(MAXIMUM_TABLE_DIRECTORY_SIZE_IN_GB);
1356  if (get_debug_level() > 0) {
1357  std::cout << "Tabular directory size is " << directory_size_in_GB << " GB\n";
1358  }
1359  if (directory_size_in_GB > 1.5 * allowed_size_in_GB) {
1360  throw DirectorySizeError(
1361  format("Maximum allowed tabular directory size is %g GB, you have exceeded 1.5 times this limit", allowed_size_in_GB));
1362  } else if (directory_size_in_GB > allowed_size_in_GB) {
1363  set_warning_string(format("Maximum allowed tabular directory size is %g GB, you have exceeded this limit", allowed_size_in_GB));
1364  }
1366  dataset->build_tables(this->AS);
1367  pack_matrices();
1368  write_tables();
1370  load_tables();
1371  // Set the flag saying tables have been successfully loaded
1372  tables_loaded = true;
1373  }
1374  }
1375  };
1376 };
1377 
1378 } /* namespace CoolProp*/
1379 
1380 #endif