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