CoolProp  6.6.0
An open-source fluid property and humid air property database
TTSEBackend.cpp
Go to the documentation of this file.
1 #if !defined(NO_TABULAR_BACKENDS)
2 
3 # include "TTSEBackend.h"
4 # include "CoolProp.h"
5 
13 double CoolProp::TTSEBackend::evaluate_single_phase_transport(SinglePhaseGriddedTableData& table, parameters output, double x, double y,
14  std::size_t i, std::size_t j) {
15  bool in_bounds = (i < table.xvec.size() - 1 && j < table.yvec.size() - 1);
16  if (!in_bounds) {
17  throw ValueError("Cell to TTSEBackend::evaluate_single_phase_transport is not valid");
18  }
19  bool is_valid = (ValidNumber(table.smolar[i][j]) && ValidNumber(table.smolar[i + 1][j]) && ValidNumber(table.smolar[i][j + 1])
20  && ValidNumber(table.smolar[i + 1][j + 1]));
21  if (!is_valid) {
22  throw ValueError("Cell to TTSEBackend::evaluate_single_phase_transport must have four valid corners for now");
23  }
24  const std::vector<std::vector<double>>& f = table.get(output);
25 
26  double x1 = table.xvec[i], x2 = table.xvec[i + 1], y1 = table.yvec[j], y2 = table.yvec[j + 1];
27  double f11 = f[i][j], f12 = f[i][j + 1], f21 = f[i + 1][j], f22 = f[i + 1][j + 1];
28  double val =
29  1 / ((x2 - x1) * (y2 - y1)) * (f11 * (x2 - x) * (y2 - y) + f21 * (x - x1) * (y2 - y) + f12 * (x2 - x) * (y - y1) + f22 * (x - x1) * (y - y1));
30 
31  // Cache the output value calculated
32  switch (output) {
33  case iconductivity:
34  _conductivity = val;
35  break;
36  case iviscosity:
37  _viscosity = val;
38  break;
39  default:
40  throw ValueError();
41  }
42  return val;
43 }
45 void CoolProp::TTSEBackend::invert_single_phase_x(const SinglePhaseGriddedTableData& table, const std::vector<std::vector<CellCoeffs>>& coeffs,
46  parameters output, double x, double y, std::size_t i, std::size_t j) {
47  connect_pointers(output, table);
48 
49  // Distances from the node
50  double deltay = y - table.yvec[j];
51 
52  // Calculate the output value desired
53  double a = 0.5 * (*d2zdx2)[i][j]; // Term multiplying deltax**2
54  double b = (*dzdx)[i][j] + deltay * (*d2zdxdy)[i][j]; // Term multiplying deltax
55  double c = (*z)[i][j] - x + deltay * (*dzdy)[i][j] + 0.5 * deltay * deltay * (*d2zdy2)[i][j];
56 
57  double deltax1 = (-b + sqrt(b * b - 4 * a * c)) / (2 * a);
58  double deltax2 = (-b - sqrt(b * b - 4 * a * c)) / (2 * a);
59 
60  // If only one is less than a multiple of x spacing, that's your solution
61  double xspacing, xratio, val;
62  if (!table.logx) {
63  xspacing = table.xvec[1] - table.xvec[0];
64  if (std::abs(deltax1) < xspacing && !(std::abs(deltax2) < xspacing)) {
65  val = deltax1 + table.xvec[i];
66  } else if (std::abs(deltax2) < xspacing && !(std::abs(deltax1) < xspacing)) {
67  val = deltax2 + table.xvec[i];
68  } else if (std::abs(deltax1) < std::abs(deltax2) && std::abs(deltax1) < 10 * xspacing) {
69  val = deltax1 + table.xvec[i];
70  } else {
71  throw ValueError(format("Cannot find the x solution; xspacing: %g dx1: %g dx2: %g", xspacing, deltax1, deltax2));
72  }
73  } else {
74  xratio = table.xvec[1] / table.xvec[0];
75  double xj = table.xvec[j];
76  double xratio1 = (xj + deltax1) / xj;
77  double xratio2 = (xj + deltax2) / xj;
78  if (xratio1 < xratio && xratio1 > 1 / xratio) {
79  val = deltax1 + table.xvec[i];
80  } else if (xratio2 < xratio && xratio2 > 1 / xratio) {
81  val = deltax2 + table.xvec[i];
82  } else if (xratio1 < xratio * 5 && xratio1 > 1 / xratio / 5) {
83  val = deltax1 + table.xvec[i];
84  } else {
85  throw ValueError(format("Cannot find the x solution; xj: %g xratio: %g xratio1: %g xratio2: %g a: %g b^2-4*a*c %g", xj, xratio, xratio1,
86  xratio2, a, b * b - 4 * a * c));
87  }
88  }
89 
90  // Cache the output value calculated
91  switch (table.xkey) {
92  case iHmolar:
93  _hmolar = val;
94  break;
95  case iT:
96  _T = val;
97  break;
98  default:
99  throw ValueError();
100  }
101 }
103 void CoolProp::TTSEBackend::invert_single_phase_y(const SinglePhaseGriddedTableData& table, const std::vector<std::vector<CellCoeffs>>& coeffs,
104  parameters output, double y, double x, std::size_t i, std::size_t j) {
105  connect_pointers(output, table);
106 
107  // Distances from the node
108  double deltax = x - table.xvec[i];
109 
110  // Calculate the output value desired
111  double a = 0.5 * (*d2zdy2)[i][j]; // Term multiplying deltay**2
112  double b = (*dzdy)[i][j] + deltax * (*d2zdxdy)[i][j]; // Term multiplying deltay
113  double c = (*z)[i][j] - y + deltax * (*dzdx)[i][j] + 0.5 * deltax * deltax * (*d2zdx2)[i][j];
114 
115  double deltay1 = (-b + sqrt(b * b - 4 * a * c)) / (2 * a);
116  double deltay2 = (-b - sqrt(b * b - 4 * a * c)) / (2 * a);
117 
118  // If only one is less than a multiple of x spacing, that's your solution
119  double yspacing, yratio, val;
120  if (!table.logy) {
121  yspacing = table.yvec[1] - table.yvec[0];
122  if (std::abs(deltay1) < yspacing && !(std::abs(deltay2) < yspacing)) {
123  val = deltay1 + table.yvec[j];
124  } else if (std::abs(deltay2) < yspacing && !(std::abs(deltay1) < yspacing)) {
125  val = deltay2 + table.yvec[j];
126  } else if (std::abs(deltay1) < std::abs(deltay2) && std::abs(deltay1) < 10 * yspacing) {
127  val = deltay1 + table.yvec[j];
128  } else {
129  throw ValueError(format("Cannot find the y solution; yspacing: %g dy1: %g dy2: %g", yspacing, deltay1, deltay2));
130  }
131  } else {
132  yratio = table.yvec[1] / table.yvec[0];
133  double yj = table.yvec[j];
134  double yratio1 = (yj + deltay1) / yj;
135  double yratio2 = (yj + deltay2) / yj;
136  if (yratio1 < yratio && yratio1 > 1 / yratio) {
137  val = deltay1 + table.yvec[j];
138  } else if (yratio2 < yratio && yratio2 > 1 / yratio) {
139  val = deltay2 + table.yvec[j];
140  } else if (std::abs(yratio1 - 1) < std::abs(yratio2 - 1)) {
141  val = deltay1 + table.yvec[j];
142  } else if (std::abs(yratio2 - 1) < std::abs(yratio1 - 1)) {
143  val = deltay2 + table.yvec[j];
144  } else {
145  throw ValueError(format("Cannot find the y solution; yj: %g yratio: %g yratio1: %g yratio2: %g a: %g b: %g b^2-4ac: %g %d %d", yj, yratio,
146  yratio1, yratio2, a, b, b * b - 4 * a * c, i, j));
147  }
148  }
149 
150  // Cache the output value calculated
151  switch (table.ykey) {
152  case iHmolar:
153  _hmolar = val;
154  break;
155  case iT:
156  _T = val;
157  break;
158  case iP:
159  _p = val;
160  break;
161  default:
162  throw ValueError();
163  }
164 }
166 double CoolProp::TTSEBackend::evaluate_single_phase(SinglePhaseGriddedTableData& table, parameters output, double x, double y, std::size_t i,
167  std::size_t j) {
168  connect_pointers(output, table);
169 
170  // Distances from the node
171  double deltax = x - table.xvec[i];
172  double deltay = y - table.yvec[j];
173 
174  // Calculate the output value desired
175  double val = (*z)[i][j] + deltax * (*dzdx)[i][j] + deltay * (*dzdy)[i][j] + 0.5 * deltax * deltax * (*d2zdx2)[i][j]
176  + 0.5 * deltay * deltay * (*d2zdy2)[i][j] + deltay * deltax * (*d2zdxdy)[i][j];
177 
178  // Cache the output value calculated
179  switch (output) {
180  case iT:
181  _T = val;
182  break;
183  case iDmolar:
184  _rhomolar = val;
185  break;
186  case iSmolar:
187  _smolar = val;
188  break;
189  case iHmolar:
190  _hmolar = val;
191  break;
192  case iUmolar:
193  _umolar = val;
194  break;
195  default:
196  throw ValueError();
197  }
198  return val;
199 }
201 double CoolProp::TTSEBackend::evaluate_single_phase_derivative(SinglePhaseGriddedTableData& table, parameters output, double x, double y,
202  std::size_t i, std::size_t j, std::size_t Nx, std::size_t Ny) {
203  if (Nx == 1 && Ny == 0) {
204  if (output == table.xkey) {
205  return 1.0;
206  }
207  if (output == table.ykey) {
208  return 0.0;
209  }
210  } else if (Ny == 1 && Nx == 0) {
211  if (output == table.ykey) {
212  return 1.0;
213  }
214  if (output == table.xkey) {
215  return 0.0;
216  }
217  }
218 
219  connect_pointers(output, table);
220 
221  // Distances from the node
222  double deltax = x - table.xvec[i];
223  double deltay = y - table.yvec[j];
224  double val;
225  // Calculate the output value desired
226  if (Nx == 1 && Ny == 0) {
227  if (output == table.xkey) {
228  return 1.0;
229  }
230  if (output == table.ykey) {
231  return 0.0;
232  }
233  val = (*dzdx)[i][j] + deltax * (*d2zdx2)[i][j] + deltay * (*d2zdxdy)[i][j];
234  } else if (Ny == 1 && Nx == 0) {
235  if (output == table.ykey) {
236  return 1.0;
237  }
238  if (output == table.xkey) {
239  return 0.0;
240  }
241  val = (*dzdy)[i][j] + deltay * (*d2zdy2)[i][j] + deltax * (*d2zdxdy)[i][j];
242  } else {
243  throw NotImplementedError("only first derivatives currently supported");
244  }
245  return val;
246 }
247 
248 #endif // !defined(NO_TABULAR_BACKENDS)
Generated on Wed Nov 29 2023 14:26:24 for CoolProp by doxygen 1.9.1