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