CoolProp  6.6.0
An open-source fluid property and humid air property database
BicubicBackend.cpp
Go to the documentation of this file.
1 #if !defined(NO_TABULAR_BACKENDS)
2 
3 # include "BicubicBackend.h"
4 # include "MatrixMath.h"
5 # include "DataStructures.h"
6 # include "Backends/Helmholtz/PhaseEnvelopeRoutines.h"
7 
8 void CoolProp::BicubicBackend::find_native_nearest_good_indices(SinglePhaseGriddedTableData& table,
9  const std::vector<std::vector<CellCoeffs>>& coeffs, double x, double y,
10  std::size_t& i, std::size_t& j) {
11  table.find_native_nearest_good_cell(x, y, i, j);
12  const CellCoeffs& cell = coeffs[i][j];
13  if (!cell.valid()) {
14  if (cell.has_valid_neighbor()) {
15  // Get new good neighbor
16  cell.get_alternate(i, j);
17  } else {
18  if (!cell.valid()) {
19  throw ValueError(format("Cell is invalid and has no good neighbors for x = %g, y= %g", x, y));
20  }
21  }
22  }
23 }
24 
26 void CoolProp::BicubicBackend::find_nearest_neighbor(SinglePhaseGriddedTableData& table, const std::vector<std::vector<CellCoeffs>>& coeffs,
27  const parameters variable1, const double value1, const parameters otherkey,
28  const double otherval, std::size_t& i, std::size_t& j) {
29  table.find_nearest_neighbor(variable1, value1, otherkey, otherval, i, j);
30  const CellCoeffs& cell = coeffs[i][j];
31  if (!cell.valid()) {
32  if (cell.has_valid_neighbor()) {
33  // Get new good neighbor
34  cell.get_alternate(i, j);
35  } else {
36  if (!cell.valid()) {
37  throw ValueError(format("Cell is invalid and has no good neighbors for x = %g, y = %g", value1, otherval));
38  }
39  }
40  }
41 }
42 
50 double CoolProp::BicubicBackend::evaluate_single_phase_transport(SinglePhaseGriddedTableData& table, parameters output, double x, double y,
51  std::size_t i, std::size_t j) {
52  // By definition i,i+1,j,j+1 are all in range and valid
53  std::vector<std::vector<double>>* f = NULL;
54  switch (output) {
55  case iconductivity:
56  f = &table.cond;
57  break;
58  case iviscosity:
59  f = &table.visc;
60  break;
61  default:
62  throw ValueError(format("invalid output variable to BicubicBackend::evaluate_single_phase_transport"));
63  }
64  double x1 = table.xvec[i], x2 = table.xvec[i + 1], y1 = table.yvec[j], y2 = table.yvec[j + 1];
65  double f11 = (*f)[i][j], f12 = (*f)[i][j + 1], f21 = (*f)[i + 1][j], f22 = (*f)[i + 1][j + 1];
66  double val =
67  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));
68 
69  // Cache the output value calculated
70  switch (output) {
71  case iconductivity:
72  _conductivity = val;
73  break;
74  case iviscosity:
75  _viscosity = val;
76  break;
77  default:
78  throw ValueError("Invalid output variable in evaluate_single_phase_transport");
79  }
80  return val;
81 }
82 // Use the single_phase table to evaluate an output
83 double CoolProp::BicubicBackend::evaluate_single_phase(const SinglePhaseGriddedTableData& table, const std::vector<std::vector<CellCoeffs>>& coeffs,
84  const parameters output, const double x, const double y, const std::size_t i,
85  const std::size_t j) {
86  // Get the cell
87  const CellCoeffs& cell = coeffs[i][j];
88 
89  // Get the alpha coefficients
90  const std::vector<double>& alpha = cell.get(output);
91 
92  // Normalized value in the range (0, 1)
93  double xhat = (x - table.xvec[i]) / (table.xvec[i + 1] - table.xvec[i]);
94  double yhat = (y - table.yvec[j]) / (table.yvec[j + 1] - table.yvec[j]);
95 
96  // Calculate the output value desired
97  // Term multiplying x^0 using Horner's method
98  double B0 = ((((0) + alpha[3 * 4 + 0]) * yhat + alpha[2 * 4 + 0]) * yhat + alpha[1 * 4 + 0]) * yhat + alpha[0 * 4 + 0];
99  // Term multiplying x^1 using Horner's method
100  double B1 = ((((0) + alpha[3 * 4 + 1]) * yhat + alpha[2 * 4 + 1]) * yhat + alpha[1 * 4 + 1]) * yhat + alpha[0 * 4 + 1];
101  // Term multiplying x^2 using Horner's method
102  double B2 = ((((0) + alpha[3 * 4 + 2]) * yhat + alpha[2 * 4 + 2]) * yhat + alpha[1 * 4 + 2]) * yhat + alpha[0 * 4 + 2];
103  // Term multiplying x^3 using Horner's method
104  double B3 = ((((0) + alpha[3 * 4 + 3]) * yhat + alpha[2 * 4 + 3]) * yhat + alpha[1 * 4 + 3]) * yhat + alpha[0 * 4 + 3];
105 
106  double val = ((((0) + B3) * xhat + B2) * xhat + B1) * xhat + B0;
107 
108  // Cache the output value calculated
109  switch (output) {
110  case iT:
111  _T = val;
112  break;
113  case iDmolar:
114  _rhomolar = val;
115  break;
116  case iSmolar:
117  _smolar = val;
118  break;
119  case iHmolar:
120  _hmolar = val;
121  break;
122  case iUmolar:
123  _umolar = val;
124  break;
125  default:
126  throw ValueError("Invalid output variable in evaluate_single_phase");
127  }
128  return val;
129 }
131 double CoolProp::BicubicBackend::evaluate_single_phase_derivative(SinglePhaseGriddedTableData& table, std::vector<std::vector<CellCoeffs>>& coeffs,
132  parameters output, double x, double y, std::size_t i, std::size_t j, std::size_t Nx,
133  std::size_t Ny) {
134 
135  // Get the cell
136  CellCoeffs& cell = coeffs[i][j];
137 
138  // Get the alpha coefficients
139  const std::vector<double>& alpha = cell.get(output);
140 
141  // Normalized value in the range (0, 1)
142  double xhat = (x - table.xvec[i]) / (table.xvec[i + 1] - table.xvec[i]);
143  double yhat = (y - table.yvec[j]) / (table.yvec[j + 1] - table.yvec[j]);
144  double dxhatdx = 1 / (table.xvec[i + 1] - table.xvec[i]);
145  double dyhatdy = 1 / (table.yvec[j + 1] - table.yvec[j]);
146 
147  // Calculate the output value desired
148  double val = 0;
149  if (Nx == 1 && Ny == 0) {
150  if (output == table.xkey) {
151  return 1.0;
152  }
153  if (output == table.ykey) {
154  return 0.0;
155  }
156  for (std::size_t l = 1; l < 4; ++l) {
157  for (std::size_t m = 0; m < 4; ++m) {
158  val += alpha[m * 4 + l] * l * pow(xhat, static_cast<int>(l - 1)) * pow(yhat, static_cast<int>(m));
159  }
160  }
161  // val is now dz/dxhat|yhat
162  return val * dxhatdx;
163  } else if (Ny == 1 && Nx == 0) {
164  if (output == table.ykey) {
165  return 1.0;
166  }
167  if (output == table.xkey) {
168  return 0.0;
169  }
170  for (std::size_t l = 0; l < 4; ++l) {
171  for (std::size_t m = 1; m < 4; ++m) {
172  val += alpha[m * 4 + l] * pow(xhat, static_cast<int>(l)) * m * pow(yhat, static_cast<int>(m - 1));
173  }
174  }
175  // val is now dz/dyhat|xhat
176  return val * dyhatdy;
177  } else {
178  throw ValueError("Invalid input");
179  }
180 }
181 
183 void CoolProp::BicubicBackend::invert_single_phase_x(const SinglePhaseGriddedTableData& table, const std::vector<std::vector<CellCoeffs>>& coeffs,
184  parameters other_key, double other, double y, std::size_t i, std::size_t j) {
185  // Get the cell
186  const CellCoeffs& cell = coeffs[i][j];
187 
188  // Get the alpha coefficients
189  const std::vector<double>& alpha = cell.get(other_key);
190 
191  // Normalized value in the range (0, 1)
192  double yhat = (y - table.yvec[j]) / (table.yvec[j + 1] - table.yvec[j]);
193 
194  double y_0 = 1, y_1 = yhat, y_2 = yhat * yhat, y_3 = yhat * yhat * yhat;
195 
196  double a = alpha[3 + 0 * 4] * y_0 + alpha[3 + 1 * 4] * y_1 + alpha[3 + 2 * 4] * y_2 + alpha[3 + 3 * 4] * y_3; // factors of xhat^3
197  double b = alpha[2 + 0 * 4] * y_0 + alpha[2 + 1 * 4] * y_1 + alpha[2 + 2 * 4] * y_2 + alpha[2 + 3 * 4] * y_3; // factors of xhar^2
198  double c = alpha[1 + 0 * 4] * y_0 + alpha[1 + 1 * 4] * y_1 + alpha[1 + 2 * 4] * y_2 + alpha[1 + 3 * 4] * y_3; // factors of xhat
199  double d = alpha[0 + 0 * 4] * y_0 + alpha[0 + 1 * 4] * y_1 + alpha[0 + 2 * 4] * y_2 + alpha[0 + 3 * 4] * y_3 - other; // constant factors
200  int N = 0;
201  double xhat0, xhat1, xhat2, val, xhat = _HUGE;
202  solve_cubic(a, b, c, d, N, xhat0, xhat1, xhat2);
203  if (N == 1) {
204  xhat = xhat0;
205  } else if (N == 2) {
206  xhat = std::abs(xhat0) < std::abs(xhat1) ? xhat0 : xhat1;
207  } else if (N == 3) {
208  if (std::abs(xhat0) < std::abs(xhat1) && std::abs(xhat0) < std::abs(xhat2)) {
209  xhat = xhat0;
210  }
211  // Already know that xhat1 < xhat0 (xhat0 is not the minimum)
212  else if (std::abs(xhat1) < std::abs(xhat2)) {
213  xhat = xhat1;
214  } else {
215  xhat = xhat2;
216  }
217  } else if (N == 0) {
218  throw ValueError("Could not find a solution in invert_single_phase_x");
219  }
220 
221  // Unpack xhat into actual value
222  // xhat = (x-x_{i})/(x_{i+1}-x_{i})
223  val = xhat * (table.xvec[i + 1] - table.xvec[i]) + table.xvec[i];
224 
225  // Cache the output value calculated
226  switch (table.xkey) {
227  case iHmolar:
228  _hmolar = val;
229  break;
230  case iT:
231  _T = val;
232  break;
233  default:
234  throw ValueError("Invalid output variable in invert_single_phase_x");
235  }
236 }
237 
239 void CoolProp::BicubicBackend::invert_single_phase_y(const SinglePhaseGriddedTableData& table, const std::vector<std::vector<CellCoeffs>>& coeffs,
240  parameters other_key, double other, double x, std::size_t i, std::size_t j) {
241  // Get the cell
242  const CellCoeffs& cell = coeffs[i][j];
243 
244  // Get the alpha coefficients
245  const std::vector<double>& alpha = cell.get(other_key);
246 
247  // Normalized value in the range (0, 1)
248  double xhat = (x - table.xvec[i]) / (table.xvec[i + 1] - table.xvec[i]);
249 
250  double x_0 = 1, x_1 = xhat, x_2 = xhat * xhat, x_3 = xhat * xhat * xhat;
251 
252  double a = alpha[0 + 3 * 4] * x_0 + alpha[1 + 3 * 4] * x_1 + alpha[2 + 3 * 4] * x_2 + alpha[3 + 3 * 4] * x_3; // factors of yhat^3 (m= 3)
253  double b = alpha[0 + 2 * 4] * x_0 + alpha[1 + 2 * 4] * x_1 + alpha[2 + 2 * 4] * x_2 + alpha[3 + 2 * 4] * x_3; // factors of yhat^2
254  double c = alpha[0 + 1 * 4] * x_0 + alpha[1 + 1 * 4] * x_1 + alpha[2 + 1 * 4] * x_2 + alpha[3 + 1 * 4] * x_3; // factors of yhat
255  double d = alpha[0 + 0 * 4] * x_0 + alpha[1 + 0 * 4] * x_1 + alpha[2 + 0 * 4] * x_2 + alpha[3 + 0 * 4] * x_3 - other; // constant factors
256  int N = 0;
257  double yhat0, yhat1, yhat2, val, yhat = _HUGE;
258  solve_cubic(a, b, c, d, N, yhat0, yhat1, yhat2);
259  if (N == 1) {
260  yhat = yhat0;
261  } else if (N == 2) {
262  yhat = std::abs(yhat0) < std::abs(yhat1) ? yhat0 : yhat1;
263  } else if (N == 3) {
264  if (std::abs(yhat0) < std::abs(yhat1) && std::abs(yhat0) < std::abs(yhat2)) {
265  yhat = yhat0;
266  }
267  // Already know that yhat1 < yhat0 (yhat0 is not the minimum)
268  else if (std::abs(yhat1) < std::abs(yhat2)) {
269  yhat = yhat1;
270  } else {
271  yhat = yhat2;
272  }
273  } else if (N == 0) {
274  throw ValueError("Could not find a solution in invert_single_phase_x");
275  }
276 
277  // Unpack xhat into actual value
278  // yhat = (y-y_{j})/(y_{j+1}-y_{j})
279  val = yhat * (table.yvec[j + 1] - table.yvec[j]) + table.yvec[j];
280 
281  // Cache the output value calculated
282  switch (table.ykey) {
283  case iP:
284  _p = val;
285  break;
286  default:
287  throw ValueError("Invalid output variable in invert_single_phase_x");
288  }
289 }
290 
291 #endif // !defined(NO_TABULAR_BACKENDS)
Generated on Wed Nov 29 2023 14:26:24 for CoolProp by doxygen 1.9.1