CoolProp 6.8.1dev
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
8void 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
26void 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
50double 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
83double 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}
131double 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
183void 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
239void 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 Jun 4 2025 04:12:26 for CoolProp by doxygen 1.9.4