CoolProp 6.8.1dev
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
13double 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}
45void 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}
103void 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}
166double 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}
201double 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 Thu Jun 5 2025 04:10:24 for CoolProp by doxygen 1.9.4