CoolProp 6.8.1dev
An open-source fluid property and humid air property database
TabularBackends.cpp
Go to the documentation of this file.
1
2#if !defined(NO_TABULAR_BACKENDS)
3
4# include "TabularBackends.h"
5# include "CoolProp.h"
6# include <sstream>
7# include "time.h"
8# include "miniz.h"
9# include <fstream>
10
13static const double Ainv_data[16 * 16] = {
14 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, -2,
15 -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, -2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
16 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, -2, -1, 0,
17 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, -2, 0, 0, 1, 1, 0, 0, -3, 0, 3, 0, 0, 0, 0, 0, -2, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
18 -3, 0, 3, 0, 0, 0, 0, 0, -2, 0, -1, 0, 9, -9, -9, 9, 6, 3, -6, -3, 6, -6, 3, -3, 4, 2, 2, 1, -6, 6, 6, -6, -3, -3, 3, 3, -4,
19 4, -2, 2, -2, -2, -1, -1, 2, 0, -2, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, -2, 0, 0, 0, 0, 0, 1, 0,
20 1, 0, -6, 6, 6, -6, -4, -2, 4, 2, -3, 3, -3, 3, -2, -1, -2, -1, 4, -4, -4, 4, 2, 2, -2, -2, 2, -2, 2, -2, 1, 1, 1, 1};
21static Eigen::Matrix<double, 16, 16> Ainv(Ainv_data);
22
23static CoolProp::TabularDataLibrary library;
24
25namespace CoolProp {
26
33template <typename T>
34void load_table(T& table, const std::string& path_to_tables, const std::string& filename) {
35
36 double tic = clock();
37 std::string path_to_table = path_to_tables + "/" + filename;
38 if (get_debug_level() > 0) {
39 std::cout << format("Loading table: %s", path_to_table.c_str()) << std::endl;
40 }
41 std::vector<char> raw;
42 try {
43 raw = get_binary_file_contents(path_to_table.c_str());
44 } catch (...) {
45 std::string err = format("Unable to load file %s", path_to_table.c_str());
46 if (get_debug_level() > 0) {
47 std::cout << "err:" << err << std::endl;
48 }
49 throw UnableToLoadError(err);
50 }
51 std::vector<unsigned char> newBuffer(raw.size() * 5);
52 uLong newBufferSize = static_cast<uLong>(newBuffer.size());
53 mz_ulong rawBufferSize = static_cast<mz_ulong>(raw.size());
54 int code;
55 do {
56 code = uncompress((unsigned char*)(&(newBuffer[0])), &newBufferSize, (unsigned char*)(&(raw[0])), rawBufferSize);
57 if (code == Z_BUF_ERROR) {
58 // Output buffer is too small, make it bigger and try again
59 newBuffer.resize(newBuffer.size() * 5);
60 newBufferSize = static_cast<uLong>(newBuffer.size());
61 } else if (code != 0) { // Something else, a big problem
62 std::string err = format("Unable to uncompress file %s with miniz code %d", path_to_table.c_str(), code);
63 if (get_debug_level() > 0) {
64 std::cout << "uncompress err:" << err << std::endl;
65 }
66 throw UnableToLoadError(err);
67 }
68 } while (code != 0);
69 // Copy the buffer from unsigned char to char (yuck)
70 std::vector<char> charbuffer(newBuffer.begin(), newBuffer.begin() + newBufferSize);
71 try {
72 msgpack::unpacked msg;
73 msgpack::unpack(msg, &(charbuffer[0]), charbuffer.size());
74 msgpack::object deserialized = msg.get();
75
76 // Call the class' deserialize function; if it is an invalid table, it will cause an exception to be thrown
77 table.deserialize(deserialized);
78 double toc = clock();
79 if (get_debug_level() > 0) {
80 std::cout << format("Loaded table: %s in %g sec.", path_to_table.c_str(), (toc - tic) / CLOCKS_PER_SEC) << std::endl;
81 }
82 } catch (std::exception& e) {
83 std::string err = format("Unable to msgpack deserialize %s; err: %s", path_to_table.c_str(), e.what());
84 if (get_debug_level() > 0) {
85 std::cout << "err: " << err << std::endl;
86 }
87 throw UnableToLoadError(err);
88 }
89}
90template <typename T>
91void write_table(const T& table, const std::string& path_to_tables, const std::string& name) {
92 msgpack::sbuffer sbuf;
93 msgpack::pack(sbuf, table);
94 std::string tabPath = std::string(path_to_tables + "/" + name + ".bin");
95 std::string zPath = tabPath + ".z";
96 std::vector<char> buffer(sbuf.size());
97 uLong outSize = static_cast<uLong>(buffer.size());
98 compress((unsigned char*)(&(buffer[0])), &outSize, (unsigned char*)(sbuf.data()), static_cast<mz_ulong>(sbuf.size()));
99 std::ofstream ofs2(zPath.c_str(), std::ofstream::binary);
100 ofs2.write(&buffer[0], outSize);
101 ofs2.close();
102
103 if (CoolProp::get_config_bool(SAVE_RAW_TABLES)) {
104 std::ofstream ofs(tabPath.c_str(), std::ofstream::binary);
105 ofs.write(sbuf.data(), sbuf.size());
106 }
107}
108
109} // namespace CoolProp
110
111void CoolProp::PureFluidSaturationTableData::build(shared_ptr<CoolProp::AbstractState>& AS) {
112 const bool debug = get_debug_level() > 5 || false;
113 if (debug) {
114 std::cout << format("***********************************************\n");
115 std::cout << format(" Saturation Table (%s) \n", AS->name().c_str());
116 std::cout << format("***********************************************\n");
117 }
118 resize(N);
119 // ------------------------
120 // Actually build the table
121 // ------------------------
122 CoolPropDbl Tmin = std::max(AS->Ttriple(), AS->Tmin());
123 AS->update(QT_INPUTS, 0, Tmin);
124 CoolPropDbl p_triple = AS->p();
125 CoolPropDbl p, pmin = p_triple, pmax = 0.9999 * AS->p_critical();
126 for (std::size_t i = 0; i < N - 1; ++i) {
127 if (i == 0) {
128 CoolProp::set_config_bool(DONT_CHECK_PROPERTY_LIMITS, true);
129 }
130 // Log spaced
131 p = exp(log(pmin) + (log(pmax) - log(pmin)) / (N - 1) * i);
132
133 // Saturated liquid
134 try {
135 AS->update(PQ_INPUTS, p, 0);
136 pL[i] = p;
137 TL[i] = AS->T();
138 rhomolarL[i] = AS->rhomolar();
139 hmolarL[i] = AS->hmolar();
140 smolarL[i] = AS->smolar();
141 umolarL[i] = AS->umolar();
142 logpL[i] = log(p);
143 logrhomolarL[i] = log(rhomolarL[i]);
144 cpmolarL[i] = AS->cpmolar();
145 cvmolarL[i] = AS->cvmolar();
146 speed_soundL[i] = AS->speed_sound();
147 } catch (std::exception& e) {
148 // That failed for some reason, go to the next pair
149 if (debug) {
150 std::cout << " " << e.what() << std::endl;
151 }
152 continue;
153 }
154 // Transport properties - if no transport properties, just keep going
155 try {
156 viscL[i] = AS->viscosity();
157 condL[i] = AS->conductivity();
158 logviscL[i] = log(viscL[i]);
159 } catch (std::exception& e) {
160 if (debug) {
161 std::cout << " " << e.what() << std::endl;
162 }
163 }
164 // Saturated vapor
165 try {
166 AS->update(PQ_INPUTS, p, 1);
167 pV[i] = p;
168 TV[i] = AS->T();
169 rhomolarV[i] = AS->rhomolar();
170 hmolarV[i] = AS->hmolar();
171 smolarV[i] = AS->smolar();
172 umolarV[i] = AS->umolar();
173 logpV[i] = log(p);
174 logrhomolarV[i] = log(rhomolarV[i]);
175 cpmolarV[i] = AS->cpmolar();
176 cvmolarV[i] = AS->cvmolar();
177 speed_soundV[i] = AS->speed_sound();
178 } catch (std::exception& e) {
179 // That failed for some reason, go to the next pair
180 if (debug) {
181 std::cout << " " << e.what() << std::endl;
182 }
183 continue;
184 }
185 // Transport properties - if no transport properties, just keep going
186 try {
187 viscV[i] = AS->viscosity();
188 condV[i] = AS->conductivity();
189 logviscV[i] = log(viscV[i]);
190 } catch (std::exception& e) {
191 if (debug) {
192 std::cout << " " << e.what() << std::endl;
193 }
194 }
195 if (i == 0) {
196 CoolProp::set_config_bool(DONT_CHECK_PROPERTY_LIMITS, false);
197 }
198 }
199 // Last point is at the critical point
200 AS->update(PQ_INPUTS, AS->p_critical(), 1);
201 std::size_t i = N - 1;
202
203 pV[i] = AS->p();
204 TV[i] = AS->T();
205 rhomolarV[i] = AS->rhomolar();
206 hmolarV[i] = AS->hmolar();
207 smolarV[i] = AS->smolar();
208 umolarV[i] = AS->umolar();
209
210 pL[i] = AS->p();
211 TL[i] = AS->T();
212 rhomolarL[i] = AS->rhomolar();
213 hmolarL[i] = AS->hmolar();
214 smolarL[i] = AS->smolar();
215 umolarL[i] = AS->umolar();
216
217 logpV[i] = log(AS->p());
218 logrhomolarV[i] = log(rhomolarV[i]);
219
220 logpL[i] = log(AS->p());
221 logrhomolarL[i] = log(rhomolarL[i]);
222}
223
224void CoolProp::SinglePhaseGriddedTableData::build(shared_ptr<CoolProp::AbstractState>& AS) {
225 CoolPropDbl x, y;
226 const bool debug = get_debug_level() > 5 || false;
227
228 resize(Nx, Ny);
229
230 if (debug) {
231 std::cout << format("***********************************************\n");
232 std::cout << format(" Single-Phase Table (%s) \n", strjoin(AS->fluid_names(), "&").c_str());
233 std::cout << format("***********************************************\n");
234 }
235 // ------------------------
236 // Actually build the table
237 // ------------------------
238 for (std::size_t i = 0; i < Nx; ++i) {
239 // Calculate the x value
240 if (logx) {
241 // Log spaced
242 x = exp(log(xmin) + (log(xmax) - log(xmin)) / (Nx - 1) * i);
243 } else {
244 // Linearly spaced
245 x = xmin + (xmax - xmin) / (Nx - 1) * i;
246 }
247 xvec[i] = x;
248 for (std::size_t j = 0; j < Ny; ++j) {
249 // Calculate the x value
250 if (logy) {
251 // Log spaced
252 y = exp(log(ymin) + (log(ymax / ymin)) / (Ny - 1) * j);
253 } else {
254 // Linearly spaced
255 y = ymin + (ymax - ymin) / (Ny - 1) * j;
256 }
257 yvec[j] = y;
258
259 if (debug) {
260 std::cout << "x: " << x << " y: " << y << std::endl;
261 }
262
263 // Generate the input pair
264 CoolPropDbl v1, v2;
265 input_pairs input_pair = generate_update_pair(xkey, x, ykey, y, v1, v2);
266
267 // --------------------
268 // Update the state
269 // --------------------
270 try {
271 AS->update(input_pair, v1, v2);
272 if (!ValidNumber(AS->rhomolar())) {
273 throw ValueError("rhomolar is invalid");
274 }
275 } catch (std::exception& e) {
276 // That failed for some reason, go to the next pair
277 if (debug) {
278 std::cout << " " << e.what() << std::endl;
279 }
280 continue;
281 }
282
283 // Skip two-phase states - they will remain as _HUGE holes in the table
284 if (is_in_closed_range(0.0, 1.0, AS->Q())) {
285 if (debug) {
286 std::cout << " 2Phase" << std::endl;
287 }
288 continue;
289 };
290
291 // --------------------
292 // State variables
293 // --------------------
294 T[i][j] = AS->T();
295 p[i][j] = AS->p();
296 rhomolar[i][j] = AS->rhomolar();
297 hmolar[i][j] = AS->hmolar();
298 smolar[i][j] = AS->smolar();
299 umolar[i][j] = AS->umolar();
300
301 // -------------------------
302 // Transport properties
303 // -------------------------
304 try {
305 visc[i][j] = AS->viscosity();
306 cond[i][j] = AS->conductivity();
307 } catch (std::exception&) {
308 // Failures will remain as holes in table
309 }
310
311 // ----------------------------------------
312 // First derivatives of state variables
313 // ----------------------------------------
314 dTdx[i][j] = AS->first_partial_deriv(iT, xkey, ykey);
315 dTdy[i][j] = AS->first_partial_deriv(iT, ykey, xkey);
316 dpdx[i][j] = AS->first_partial_deriv(iP, xkey, ykey);
317 dpdy[i][j] = AS->first_partial_deriv(iP, ykey, xkey);
318 drhomolardx[i][j] = AS->first_partial_deriv(iDmolar, xkey, ykey);
319 drhomolardy[i][j] = AS->first_partial_deriv(iDmolar, ykey, xkey);
320 dhmolardx[i][j] = AS->first_partial_deriv(iHmolar, xkey, ykey);
321 dhmolardy[i][j] = AS->first_partial_deriv(iHmolar, ykey, xkey);
322 dsmolardx[i][j] = AS->first_partial_deriv(iSmolar, xkey, ykey);
323 dsmolardy[i][j] = AS->first_partial_deriv(iSmolar, ykey, xkey);
324 dumolardx[i][j] = AS->first_partial_deriv(iUmolar, xkey, ykey);
325 dumolardy[i][j] = AS->first_partial_deriv(iUmolar, ykey, xkey);
326
327 // ----------------------------------------
328 // Second derivatives of state variables
329 // ----------------------------------------
330 d2Tdx2[i][j] = AS->second_partial_deriv(iT, xkey, ykey, xkey, ykey);
331 d2Tdxdy[i][j] = AS->second_partial_deriv(iT, xkey, ykey, ykey, xkey);
332 d2Tdy2[i][j] = AS->second_partial_deriv(iT, ykey, xkey, ykey, xkey);
333 d2pdx2[i][j] = AS->second_partial_deriv(iP, xkey, ykey, xkey, ykey);
334 d2pdxdy[i][j] = AS->second_partial_deriv(iP, xkey, ykey, ykey, xkey);
335 d2pdy2[i][j] = AS->second_partial_deriv(iP, ykey, xkey, ykey, xkey);
336 d2rhomolardx2[i][j] = AS->second_partial_deriv(iDmolar, xkey, ykey, xkey, ykey);
337 d2rhomolardxdy[i][j] = AS->second_partial_deriv(iDmolar, xkey, ykey, ykey, xkey);
338 d2rhomolardy2[i][j] = AS->second_partial_deriv(iDmolar, ykey, xkey, ykey, xkey);
339 d2hmolardx2[i][j] = AS->second_partial_deriv(iHmolar, xkey, ykey, xkey, ykey);
340 d2hmolardxdy[i][j] = AS->second_partial_deriv(iHmolar, xkey, ykey, ykey, xkey);
341 d2hmolardy2[i][j] = AS->second_partial_deriv(iHmolar, ykey, xkey, ykey, xkey);
342 d2smolardx2[i][j] = AS->second_partial_deriv(iSmolar, xkey, ykey, xkey, ykey);
343 d2smolardxdy[i][j] = AS->second_partial_deriv(iSmolar, xkey, ykey, ykey, xkey);
344 d2smolardy2[i][j] = AS->second_partial_deriv(iSmolar, ykey, xkey, ykey, xkey);
345 d2umolardx2[i][j] = AS->second_partial_deriv(iUmolar, xkey, ykey, xkey, ykey);
346 d2umolardxdy[i][j] = AS->second_partial_deriv(iUmolar, xkey, ykey, ykey, xkey);
347 d2umolardy2[i][j] = AS->second_partial_deriv(iUmolar, ykey, xkey, ykey, xkey);
348 }
349 }
350}
351std::string CoolProp::TabularBackend::path_to_tables(void) {
352 std::vector<std::string> fluids = AS->fluid_names();
353 std::vector<CoolPropDbl> fractions = AS->get_mole_fractions();
354 std::vector<std::string> components;
355 for (std::size_t i = 0; i < fluids.size(); ++i) {
356 components.push_back(format("%s[%0.10Lf]", fluids[i].c_str(), fractions[i]));
357 }
358 std::string table_directory = get_home_dir() + "/.CoolProp/Tables/";
359 std::string alt_table_directory = get_config_string(ALTERNATIVE_TABLES_DIRECTORY);
360 if (!alt_table_directory.empty()) {
361 table_directory = alt_table_directory;
362 }
363 return table_directory + AS->backend_name() + "(" + strjoin(components, "&") + ")";
364}
365
366void CoolProp::TabularBackend::write_tables() {
367 std::string path_to_tables = this->path_to_tables();
368 make_dirs(path_to_tables);
369 bool loaded = false;
370 dataset = library.get_set_of_tables(this->AS, loaded);
371 PackablePhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
372 PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
373 SinglePhaseGriddedTableData& single_phase_logph = dataset->single_phase_logph;
374 SinglePhaseGriddedTableData& single_phase_logpT = dataset->single_phase_logpT;
375 write_table(single_phase_logph, path_to_tables, "single_phase_logph");
376 write_table(single_phase_logpT, path_to_tables, "single_phase_logpT");
377 write_table(pure_saturation, path_to_tables, "pure_saturation");
378 write_table(phase_envelope, path_to_tables, "phase_envelope");
379}
380void CoolProp::TabularBackend::load_tables() {
381 bool loaded = false;
382 dataset = library.get_set_of_tables(this->AS, loaded);
383 if (loaded == false) {
384 throw UnableToLoadError("Could not load tables");
385 }
386 if (get_debug_level() > 0) {
387 std::cout << "Tables loaded" << std::endl;
388 }
389}
390
391CoolPropDbl CoolProp::TabularBackend::calc_saturated_vapor_keyed_output(parameters key) {
392 PhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
393 PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
394 double factor = 1.0;
395 mass_to_molar(key, factor, molar_mass());
396 if (is_mixture) {
397 return phase_envelope_sat(phase_envelope, key, iP, _p) * factor;
398 } else {
399 return pure_saturation.evaluate(key, _p, 1, cached_saturation_iL, cached_saturation_iV) * factor;
400 }
401}
402CoolPropDbl CoolProp::TabularBackend::calc_saturated_liquid_keyed_output(parameters key) {
403 PhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
404 PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
405 double factor = 1.0;
406 mass_to_molar(key, factor, molar_mass());
407 if (is_mixture) {
408 return phase_envelope_sat(phase_envelope, key, iP, _p) * factor;
409 } else {
410 return pure_saturation.evaluate(key, _p, 0, cached_saturation_iL, cached_saturation_iV) * factor;
411 }
412};
413
414CoolPropDbl CoolProp::TabularBackend::calc_p(void) {
415 PhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
416 if (using_single_phase_table) {
417 return _p;
418 } else {
419 if (is_mixture) {
420 return phase_envelope_sat(phase_envelope, iP, iT, _T);
421 } else {
422 return _p;
423 }
424 }
425}
426CoolPropDbl CoolProp::TabularBackend::calc_T(void) {
427 PhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
428 PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
429 if (using_single_phase_table) {
430 switch (selected_table) {
431 case SELECTED_PH_TABLE:
432 return evaluate_single_phase_phmolar(iT, cached_single_phase_i, cached_single_phase_j);
433 case SELECTED_PT_TABLE:
434 return _T;
435 case SELECTED_NO_TABLE:
436 throw ValueError("table not selected");
437 }
438 return _HUGE; // not needed, will never be hit, just to make compiler happy
439 } else {
440 if (is_mixture) {
441 return phase_envelope_sat(phase_envelope, iT, iP, _p);
442 } else {
443 if (ValidNumber(_T)) {
444 return _T;
445 } else {
446 return pure_saturation.evaluate(iT, _p, _Q, cached_saturation_iL, cached_saturation_iV);
447 }
448 }
449 }
450}
451CoolPropDbl CoolProp::TabularBackend::calc_rhomolar(void) {
452 PhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
453 PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
454 if (using_single_phase_table) {
455 switch (selected_table) {
456 case SELECTED_PH_TABLE:
457 return evaluate_single_phase_phmolar(iDmolar, cached_single_phase_i, cached_single_phase_j);
458 case SELECTED_PT_TABLE:
459 return evaluate_single_phase_pT(iDmolar, cached_single_phase_i, cached_single_phase_j);
460 case SELECTED_NO_TABLE:
461 throw ValueError("table not selected");
462 }
463 return _HUGE; // not needed, will never be hit, just to make compiler happy
464 } else {
465 if (is_mixture) {
466 return phase_envelope_sat(phase_envelope, iDmolar, iP, _p);
467 } else {
468 return pure_saturation.evaluate(iDmolar, _p, _Q, cached_saturation_iL, cached_saturation_iV);
469 }
470 }
471}
472CoolPropDbl CoolProp::TabularBackend::calc_hmolar(void) {
473 PhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
474 PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
475 if (using_single_phase_table) {
476 switch (selected_table) {
477 case SELECTED_PH_TABLE:
478 return _hmolar;
479 case SELECTED_PT_TABLE:
480 return evaluate_single_phase_pT(iHmolar, cached_single_phase_i, cached_single_phase_j);
481 case SELECTED_NO_TABLE:
482 throw ValueError("table not selected");
483 }
484 return _HUGE; // not needed, will never be hit, just to make compiler happy
485 } else {
486 if (is_mixture) {
487 return phase_envelope_sat(phase_envelope, iHmolar, iP, _p);
488 } else {
489 return pure_saturation.evaluate(iHmolar, _p, _Q, cached_saturation_iL, cached_saturation_iV);
490 }
491 }
492}
493CoolPropDbl CoolProp::TabularBackend::calc_smolar(void) {
494 PhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
495 PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
496 if (using_single_phase_table) {
497 switch (selected_table) {
498 case SELECTED_PH_TABLE:
499 return evaluate_single_phase_phmolar(iSmolar, cached_single_phase_i, cached_single_phase_j);
500 case SELECTED_PT_TABLE:
501 return evaluate_single_phase_pT(iSmolar, cached_single_phase_i, cached_single_phase_j);
502 case SELECTED_NO_TABLE:
503 throw ValueError("table not selected");
504 }
505 return _HUGE; // not needed, will never be hit, just to make compiler happy
506 } else {
507 if (is_mixture) {
508 return phase_envelope_sat(phase_envelope, iSmolar, iP, _p);
509 } else {
510 return pure_saturation.evaluate(iSmolar, _p, _Q, cached_saturation_iL, cached_saturation_iV);
511 }
512 }
513}
514CoolPropDbl CoolProp::TabularBackend::calc_umolar(void) {
515 PhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
516 PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
517 if (using_single_phase_table) {
518 switch (selected_table) {
519 case SELECTED_PH_TABLE:
520 return evaluate_single_phase_phmolar(iUmolar, cached_single_phase_i, cached_single_phase_j);
521 case SELECTED_PT_TABLE:
522 return evaluate_single_phase_pT(iUmolar, cached_single_phase_i, cached_single_phase_j);
523 case SELECTED_NO_TABLE:
524 throw ValueError("table not selected");
525 }
526 return _HUGE; // not needed, will never be hit, just to make compiler happy
527 } else {
528 if (is_mixture) {
529 // h = u + pv -> u = h - pv
530 CoolPropDbl hmolar = phase_envelope_sat(phase_envelope, iHmolar, iP, _p);
531 CoolPropDbl rhomolar = phase_envelope_sat(phase_envelope, iDmolar, iP, _p);
532 return hmolar - _p / rhomolar;
533 } else {
534 return pure_saturation.evaluate(iUmolar, _p, _Q, cached_saturation_iL, cached_saturation_iV);
535 }
536 }
537}
538CoolPropDbl CoolProp::TabularBackend::calc_cpmolar(void) {
539 PhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
540 PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
541 if (using_single_phase_table) {
542 return calc_first_partial_deriv(iHmolar, iT, iP);
543 } else {
544 if (is_mixture) {
545 return phase_envelope_sat(phase_envelope, iCpmolar, iP, _p);
546 } else {
547 return pure_saturation.evaluate(iCpmolar, _p, _Q, cached_saturation_iL, cached_saturation_iV);
548 }
549 }
550}
551CoolPropDbl CoolProp::TabularBackend::calc_cvmolar(void) {
552 PhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
553 PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
554 if (using_single_phase_table) {
555 return calc_first_partial_deriv(iUmolar, iT, iDmolar);
556 } else {
557 if (is_mixture) {
558 return phase_envelope_sat(phase_envelope, iCvmolar, iP, _p);
559 } else {
560 return pure_saturation.evaluate(iCvmolar, _p, _Q, cached_saturation_iL, cached_saturation_iV);
561 }
562 }
563}
564
565CoolPropDbl CoolProp::TabularBackend::calc_viscosity(void) {
566 PhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
567 PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
568 if (using_single_phase_table) {
569 switch (selected_table) {
570 case SELECTED_PH_TABLE:
571 return evaluate_single_phase_phmolar_transport(iviscosity, cached_single_phase_i, cached_single_phase_j);
572 case SELECTED_PT_TABLE:
573 return evaluate_single_phase_pT_transport(iviscosity, cached_single_phase_i, cached_single_phase_j);
574 case SELECTED_NO_TABLE:
575 throw ValueError("table not selected");
576 }
577 return _HUGE; // not needed, will never be hit, just to make compiler happy
578 } else {
579 if (is_mixture) {
580 return phase_envelope_sat(phase_envelope, iviscosity, iP, _p);
581 } else {
582 return pure_saturation.evaluate(iviscosity, _p, _Q, cached_saturation_iL, cached_saturation_iV);
583 }
584 }
585}
586CoolPropDbl CoolProp::TabularBackend::calc_conductivity(void) {
587 PhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
588 PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
589 if (using_single_phase_table) {
590 switch (selected_table) {
591 case SELECTED_PH_TABLE:
592 return evaluate_single_phase_phmolar_transport(iconductivity, cached_single_phase_i, cached_single_phase_j);
593 case SELECTED_PT_TABLE:
594 return evaluate_single_phase_pT_transport(iconductivity, cached_single_phase_i, cached_single_phase_j);
595 case SELECTED_NO_TABLE:
596 throw ValueError("table not selected");
597 }
598 return _HUGE; // not needed, will never be hit, just to make compiler happy
599 } else {
600 if (is_mixture) {
601 return phase_envelope_sat(phase_envelope, iconductivity, iP, _p);
602 } else {
603 return pure_saturation.evaluate(iconductivity, _p, _Q, cached_saturation_iL, cached_saturation_iV);
604 }
605 }
606}
607CoolPropDbl CoolProp::TabularBackend::calc_speed_sound(void) {
608 PhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
609 PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
610 if (using_single_phase_table) {
611 return sqrt(1 / molar_mass() * cpmolar() / cvmolar() * first_partial_deriv(iP, iDmolar, iT));
612 } else {
613 if (is_mixture) {
614 return phase_envelope_sat(phase_envelope, ispeed_sound, iP, _p);
615 } else {
616 return pure_saturation.evaluate(ispeed_sound, _p, _Q, cached_saturation_iL, cached_saturation_iV);
617 }
618 }
619}
620CoolPropDbl CoolProp::TabularBackend::calc_first_partial_deriv(parameters Of, parameters Wrt, parameters Constant) {
621 if (using_single_phase_table) {
622 CoolPropDbl dOf_dx, dOf_dy, dWrt_dx, dWrt_dy, dConstant_dx, dConstant_dy;
623
624 // If a mass-based parameter is provided, get a conversion factor and change the key to the molar-based key
625 double Of_conversion_factor = 1.0, Wrt_conversion_factor = 1.0, Constant_conversion_factor = 1.0, MM = AS->molar_mass();
626 mass_to_molar(Of, Of_conversion_factor, MM);
627 mass_to_molar(Wrt, Wrt_conversion_factor, MM);
628 mass_to_molar(Constant, Constant_conversion_factor, MM);
629
630 switch (selected_table) {
631 case SELECTED_PH_TABLE: {
632 dOf_dx = evaluate_single_phase_phmolar_derivative(Of, cached_single_phase_i, cached_single_phase_j, 1, 0);
633 dOf_dy = evaluate_single_phase_phmolar_derivative(Of, cached_single_phase_i, cached_single_phase_j, 0, 1);
634 dWrt_dx = evaluate_single_phase_phmolar_derivative(Wrt, cached_single_phase_i, cached_single_phase_j, 1, 0);
635 dWrt_dy = evaluate_single_phase_phmolar_derivative(Wrt, cached_single_phase_i, cached_single_phase_j, 0, 1);
636 dConstant_dx = evaluate_single_phase_phmolar_derivative(Constant, cached_single_phase_i, cached_single_phase_j, 1, 0);
637 dConstant_dy = evaluate_single_phase_phmolar_derivative(Constant, cached_single_phase_i, cached_single_phase_j, 0, 1);
638 break;
639 }
640 case SELECTED_PT_TABLE: {
641 dOf_dx = evaluate_single_phase_pT_derivative(Of, cached_single_phase_i, cached_single_phase_j, 1, 0);
642 dOf_dy = evaluate_single_phase_pT_derivative(Of, cached_single_phase_i, cached_single_phase_j, 0, 1);
643 dWrt_dx = evaluate_single_phase_pT_derivative(Wrt, cached_single_phase_i, cached_single_phase_j, 1, 0);
644 dWrt_dy = evaluate_single_phase_pT_derivative(Wrt, cached_single_phase_i, cached_single_phase_j, 0, 1);
645 dConstant_dx = evaluate_single_phase_pT_derivative(Constant, cached_single_phase_i, cached_single_phase_j, 1, 0);
646 dConstant_dy = evaluate_single_phase_pT_derivative(Constant, cached_single_phase_i, cached_single_phase_j, 0, 1);
647 break;
648 }
649 case SELECTED_NO_TABLE:
650 throw ValueError("table not selected");
651 }
652 double val = (dOf_dx * dConstant_dy - dOf_dy * dConstant_dx) / (dWrt_dx * dConstant_dy - dWrt_dy * dConstant_dx);
653 return val * Of_conversion_factor / Wrt_conversion_factor;
654 } else {
655 throw ValueError(format("Inputs [rho: %g mol/m3, T: %g K, p: %g Pa] are two-phase; cannot use single-phase derivatives", _rhomolar, _T, _p));
656 }
657};
658
659CoolPropDbl CoolProp::TabularBackend::calc_first_saturation_deriv(parameters Of1, parameters Wrt1) {
660 PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
661 if (AS->get_mole_fractions().size() > 1) {
662 throw ValueError("calc_first_saturation_deriv not available for mixtures");
663 }
664 if (std::abs(_Q) < 1e-6) {
665 return pure_saturation.first_saturation_deriv(Of1, Wrt1, 0, keyed_output(Wrt1), cached_saturation_iL);
666 } else if (std::abs(_Q - 1) < 1e-6) {
667 return pure_saturation.first_saturation_deriv(Of1, Wrt1, 1, keyed_output(Wrt1), cached_saturation_iV);
668 } else {
669 throw ValueError(format("Quality [%Lg] must be either 0 or 1 to within 1 ppm", _Q));
670 }
671}
672CoolPropDbl CoolProp::TabularBackend::calc_first_two_phase_deriv(parameters Of, parameters Wrt, parameters Constant) {
673 PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
674 if (Of == iDmolar && Wrt == iHmolar && Constant == iP) {
675 CoolPropDbl rhoL = pure_saturation.evaluate(iDmolar, _p, 0, cached_saturation_iL, cached_saturation_iV);
676 CoolPropDbl rhoV = pure_saturation.evaluate(iDmolar, _p, 1, cached_saturation_iL, cached_saturation_iV);
677 CoolPropDbl hL = pure_saturation.evaluate(iHmolar, _p, 0, cached_saturation_iL, cached_saturation_iV);
678 CoolPropDbl hV = pure_saturation.evaluate(iHmolar, _p, 1, cached_saturation_iL, cached_saturation_iV);
679 return -POW2(rhomolar()) * (1 / rhoV - 1 / rhoL) / (hV - hL);
680 } else if (Of == iDmass && Wrt == iHmass && Constant == iP) {
681 return first_two_phase_deriv(iDmolar, iHmolar, iP) * POW2(molar_mass());
682 } else if (Of == iDmolar && Wrt == iP && Constant == iHmolar) {
683 // v = 1/rho; drhodv = -rho^2; dvdrho = -1/rho^2
684 CoolPropDbl rhoL = pure_saturation.evaluate(iDmolar, _p, 0, cached_saturation_iL, cached_saturation_iV);
685 CoolPropDbl rhoV = pure_saturation.evaluate(iDmolar, _p, 1, cached_saturation_iL, cached_saturation_iV);
686 CoolPropDbl hL = pure_saturation.evaluate(iHmolar, _p, 0, cached_saturation_iL, cached_saturation_iV);
687 CoolPropDbl hV = pure_saturation.evaluate(iHmolar, _p, 1, cached_saturation_iL, cached_saturation_iV);
688 CoolPropDbl dvdrhoL = -1 / POW2(rhoL);
689 CoolPropDbl dvdrhoV = -1 / POW2(rhoV);
690 CoolPropDbl dvL_dp = dvdrhoL * pure_saturation.first_saturation_deriv(iDmolar, iP, 0, _p, cached_saturation_iL);
691 CoolPropDbl dvV_dp = dvdrhoV * pure_saturation.first_saturation_deriv(iDmolar, iP, 1, _p, cached_saturation_iV);
692 CoolPropDbl dhL_dp = pure_saturation.first_saturation_deriv(iHmolar, iP, 0, _p, cached_saturation_iL);
693 CoolPropDbl dhV_dp = pure_saturation.first_saturation_deriv(iHmolar, iP, 1, _p, cached_saturation_iV);
694 CoolPropDbl dxdp_h = (Q() * dhV_dp + (1 - Q()) * dhL_dp) / (hL - hV);
695 CoolPropDbl dvdp_h = dvL_dp + dxdp_h * (1 / rhoV - 1 / rhoL) + Q() * (dvV_dp - dvL_dp);
696 return -POW2(rhomolar()) * dvdp_h;
697 } else if (Of == iDmass && Wrt == iP && Constant == iHmass) {
698 return first_two_phase_deriv(iDmolar, iP, iHmolar) * molar_mass();
699 } else {
700 throw ValueError("These inputs are not supported to calc_first_two_phase_deriv");
701 }
702}
703
704CoolPropDbl CoolProp::TabularBackend::calc_first_two_phase_deriv_splined(parameters Of, parameters Wrt, parameters Constant, CoolPropDbl x_end) {
705 // Note: If you need all three values (drho_dh__p, drho_dp__h and rho_spline),
706 // you should calculate drho_dp__h first to avoid duplicate calculations.
707
708 bool drho_dh__p = false;
709 bool drho_dp__h = false;
710 bool rho_spline = false;
711
712 if (Of == iDmolar && Wrt == iHmolar && Constant == iP) {
713 drho_dh__p = true;
714 if (_drho_spline_dh__constp) return _drho_spline_dh__constp;
715 } else if (Of == iDmass && Wrt == iHmass && Constant == iP) {
716 return first_two_phase_deriv_splined(iDmolar, iHmolar, iP, x_end) * POW2(molar_mass());
717 } else if (Of == iDmolar && Wrt == iP && Constant == iHmolar) {
718 drho_dp__h = true;
719 if (_drho_spline_dp__consth) return _drho_spline_dp__consth;
720 } else if (Of == iDmass && Wrt == iP && Constant == iHmass) {
721 return first_two_phase_deriv_splined(iDmolar, iP, iHmolar, x_end) * molar_mass();
722 }
723 // Add the special case for the splined density
724 else if (Of == iDmolar && Wrt == iDmolar && Constant == iDmolar) {
725 rho_spline = true;
726 if (_rho_spline) return _rho_spline;
727 } else if (Of == iDmass && Wrt == iDmass && Constant == iDmass) {
728 return first_two_phase_deriv_splined(iDmolar, iDmolar, iDmolar, x_end) * molar_mass();
729 } else {
730 throw ValueError("These inputs are not supported to calc_first_two_phase_deriv");
731 }
732
733 if (_Q > x_end) {
734 throw ValueError(format("Q [%g] is greater than x_end [%Lg]", _Q, x_end).c_str());
735 }
736 if (_phase != iphase_twophase) {
737 throw ValueError(format("state is not two-phase"));
738 }
739
740 // TODO: replace AS with a cloned instance of TTSE or BICUBIC, add "clone()" as mandatory function in base class
741 //shared_ptr<TabularBackend>
742 // Liq(this->clone())),
743 // End(this->clone()));
744 //
745 //Liq->specify_phase(iphase_liquid);
746 //Liq->_Q = -1;
747 //Liq->update_DmolarT_direct(SatL->rhomolar(), SatL->T());
748 //End->update(QT_INPUTS, x_end, SatL->T());
749
750 // Start with quantities needed for all calculations
751 PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
752 CoolPropDbl hL = pure_saturation.evaluate(iHmolar, _p, 0, cached_saturation_iL, cached_saturation_iV);
753 CoolPropDbl hV = pure_saturation.evaluate(iHmolar, _p, 1, cached_saturation_iL, cached_saturation_iV);
754 CoolPropDbl hE = pure_saturation.evaluate(iHmolar, _p, x_end, cached_saturation_iL, cached_saturation_iV);
755
756 CoolPropDbl dL = pure_saturation.evaluate(iDmolar, _p, 0, cached_saturation_iL, cached_saturation_iV);
757 CoolPropDbl dV = pure_saturation.evaluate(iDmolar, _p, 1, cached_saturation_iL, cached_saturation_iV);
758 CoolPropDbl dE = pure_saturation.evaluate(iDmolar, _p, x_end, cached_saturation_iL, cached_saturation_iV);
759
760 CoolPropDbl Delta = Q() * (hV - hL);
761 CoolPropDbl Delta_end = hE - hL;
762
763 CoolPropDbl TL = pure_saturation.evaluate(iT, _p, 0, cached_saturation_iL, cached_saturation_iV);
764 CoolPropDbl TE = pure_saturation.evaluate(iT, _p, x_end, cached_saturation_iL, cached_saturation_iV);
765
766 // TODO: We cheat here and fake access to a derived class.
767 // Liquid state
768 AS->specify_phase(iphase_liquid);
769 AS->update(DmolarT_INPUTS, dL, TL);
770 CoolPropDbl drho_dh_liq__constp = AS->first_partial_deriv(iDmolar, iHmolar, iP);
771 CoolPropDbl d2rhodhdp_liq = AS->second_partial_deriv(iDmolar, iHmolar, iP, iP, iHmolar);
772 // End of spline
773 AS->specify_phase(iphase_twophase);
774 AS->update(DmolarT_INPUTS, dE, TE);
775 CoolPropDbl drho_dh_end__constp = AS->first_partial_deriv(iDmolar, iHmolar, iP);
776 CoolPropDbl d2rhodhdp_end = AS->second_partial_deriv(iDmolar, iHmolar, iP, iP, iHmolar);
777
778 // Spline coordinates a, b, c, d
779 CoolPropDbl Abracket = (2 * dL - 2 * dE + Delta_end * (drho_dh_liq__constp + drho_dh_end__constp));
780 CoolPropDbl a = 1 / POW3(Delta_end) * Abracket;
781 CoolPropDbl b = 3 / POW2(Delta_end) * (-dL + dE) - 1 / Delta_end * (drho_dh_end__constp + 2 * drho_dh_liq__constp);
782 CoolPropDbl c = drho_dh_liq__constp;
783 CoolPropDbl d = dL;
784
785 // Either the spline value or drho/dh|p can be directly evaluated now
786 _rho_spline = a * POW3(Delta) + b * POW2(Delta) + c * Delta + d;
787 _drho_spline_dh__constp = 3 * a * POW2(Delta) + 2 * b * Delta + c;
788 if (rho_spline) return _rho_spline;
789 if (drho_dh__p) return _drho_spline_dh__constp;
790
791 // It's drho/dp|h
792 // ... calculate some more things
793
794 // Derivatives *along* the saturation curve using the special internal method
795 CoolPropDbl dhL_dp_sat = pure_saturation.first_saturation_deriv(iHmolar, iP, 0, _p, cached_saturation_iL);
796 CoolPropDbl dhV_dp_sat = pure_saturation.first_saturation_deriv(iHmolar, iP, 1, _p, cached_saturation_iV);
797 CoolPropDbl drhoL_dp_sat = pure_saturation.first_saturation_deriv(iDmolar, iP, 0, _p, cached_saturation_iL);
798 CoolPropDbl drhoV_dp_sat = pure_saturation.first_saturation_deriv(iDmolar, iP, 1, _p, cached_saturation_iV);
799 //CoolPropDbl rhoV = SatV->keyed_output(rho_key);
800 //CoolPropDbl rhoL = SatL->keyed_output(rho_key);
801 CoolPropDbl drho_dp_end = POW2(dE) * (x_end / POW2(dV) * drhoV_dp_sat + (1 - x_end) / POW2(dL) * drhoL_dp_sat);
802
803 // Faking single-phase
804 //CoolPropDbl drho_dp__consth_liq = Liq->first_partial_deriv(rho_key, p_key, h_key);
805 //CoolPropDbl d2rhodhdp_liq = Liq->second_partial_deriv(rho_key, h_key, p_key, p_key, h_key); // ?
806
807 // Derivatives at the end point
808 // CoolPropDbl drho_dp__consth_end = End->calc_first_two_phase_deriv(rho_key, p_key, h_key);
809 //CoolPropDbl d2rhodhdp_end = End->calc_second_two_phase_deriv(rho_key, h_key, p_key, p_key, h_key);
810
811 // Reminder:
812 // Delta = Q()*(hV-hL) = h-hL
813 // Delta_end = x_end*(hV-hL);
814 CoolPropDbl d_Delta_dp__consth = -dhL_dp_sat;
815 CoolPropDbl d_Delta_end_dp__consth = x_end * (dhV_dp_sat - dhL_dp_sat);
816
817 // First pressure derivative at constant h of the coefficients a,b,c,d
818 // CoolPropDbl Abracket = (2*rho_liq - 2*rho_end + Delta_end * (drho_dh_liq__constp + drho_dh_end));
819 CoolPropDbl d_Abracket_dp_consth = (2 * drhoL_dp_sat - 2 * drho_dp_end + Delta_end * (d2rhodhdp_liq + d2rhodhdp_end)
820 + d_Delta_end_dp__consth * (drho_dh_liq__constp + drho_dh_end__constp));
821 CoolPropDbl da_dp = 1 / POW3(Delta_end) * d_Abracket_dp_consth + Abracket * (-3 / POW4(Delta_end) * d_Delta_end_dp__consth);
822 CoolPropDbl db_dp = -6 / POW3(Delta_end) * d_Delta_end_dp__consth * (dE - dL) + 3 / POW2(Delta_end) * (drho_dp_end - drhoL_dp_sat)
823 + (1 / POW2(Delta_end) * d_Delta_end_dp__consth) * (drho_dh_end__constp + 2 * drho_dh_liq__constp)
824 - (1 / Delta_end) * (d2rhodhdp_end + 2 * d2rhodhdp_liq);
825 CoolPropDbl dc_dp = d2rhodhdp_liq;
826 CoolPropDbl dd_dp = drhoL_dp_sat;
827
828 _drho_spline_dp__consth =
829 (3 * a * POW2(Delta) + 2 * b * Delta + c) * d_Delta_dp__consth + POW3(Delta) * da_dp + POW2(Delta) * db_dp + Delta * dc_dp + dd_dp;
830 if (drho_dp__h) return _drho_spline_dp__consth;
831
832 throw ValueError("Something went wrong in TabularBackend::calc_first_two_phase_deriv_splined");
833 return _HUGE;
834}
835
836void CoolProp::TabularBackend::update(CoolProp::input_pairs input_pair, double val1, double val2) {
837
838 if (get_debug_level() > 0) {
839 std::cout << format("update(%s,%g,%g)\n", get_input_pair_short_desc(input_pair).c_str(), val1, val2);
840 }
841
842 // Clear cached variables
843 clear();
844
845 // Convert to mass-based units if necessary
846 CoolPropDbl ld_value1 = val1, ld_value2 = val2;
847 mass_to_molar_inputs(input_pair, ld_value1, ld_value2);
848 val1 = ld_value1;
849 val2 = ld_value2;
850
851 // Check the tables, build if necessary
852 check_tables();
853
854 // Flush the cached indices (set to large number)
855 cached_single_phase_i = std::numeric_limits<std::size_t>::max();
856 cached_single_phase_j = std::numeric_limits<std::size_t>::max();
857 cached_saturation_iL = std::numeric_limits<std::size_t>::max();
858 cached_saturation_iV = std::numeric_limits<std::size_t>::max();
859
860 // To start, set quality to value that is impossible
861 _Q = -1000;
862
863 PureFluidSaturationTableData& pure_saturation = dataset->pure_saturation;
864 PhaseEnvelopeData& phase_envelope = dataset->phase_envelope;
865 SinglePhaseGriddedTableData& single_phase_logph = dataset->single_phase_logph;
866 SinglePhaseGriddedTableData& single_phase_logpT = dataset->single_phase_logpT;
867
868 switch (input_pair) {
869 case HmolarP_INPUTS: {
870 _hmolar = val1;
871 _p = val2;
872 if (!single_phase_logph.native_inputs_are_in_range(_hmolar, _p)) {
873 // Use the AbstractState instance
874 using_single_phase_table = false;
875 if (get_debug_level() > 5) {
876 std::cout << "inputs are not in range";
877 }
878 throw ValueError(format("inputs are not in range, hmolar=%Lg, p=%Lg", static_cast<CoolPropDbl>(_hmolar), _p));
879 } else {
880 using_single_phase_table = true; // Use the table !
881 std::size_t iL = std::numeric_limits<std::size_t>::max(), iV = std::numeric_limits<std::size_t>::max(), iclosest = 0;
882 CoolPropDbl hL = 0, hV = 0;
883 SimpleState closest_state;
884 bool is_two_phase = false;
885 // If phase is imposed, use it, but only if it's single phase:
886 // - Imposed two phase still needs to determine saturation limits by calling pure_saturation.is_inside().
887 // - There's no speed increase to be gained by imposing two phase.
888 if ((imposed_phase_index == iphase_not_imposed) || (imposed_phase_index == iphase_twophase)) {
889 if (is_mixture) {
890 is_two_phase = PhaseEnvelopeRoutines::is_inside(phase_envelope, iP, _p, iHmolar, _hmolar, iclosest, closest_state);
891 } else {
892 is_two_phase = pure_saturation.is_inside(iP, _p, iHmolar, _hmolar, iL, iV, hL, hV);
893 }
894 }
895 // Phase determined or imposed, now interpolate results
896 if (is_two_phase) {
897 using_single_phase_table = false;
898 _Q = (static_cast<double>(_hmolar) - hL) / (hV - hL);
899 if (!is_in_closed_range(0.0, 1.0, static_cast<double>(_Q))) {
900 throw ValueError(
901 format("vapor quality is not in (0,1) for hmolar: %g p: %g, hL: %g hV: %g ", static_cast<double>(_hmolar), _p, hL, hV));
902 } else {
903 cached_saturation_iL = iL;
904 cached_saturation_iV = iV;
905 _phase = iphase_twophase;
906 }
907 } else {
908 selected_table = SELECTED_PH_TABLE;
909 // Find and cache the indices i, j
910 find_native_nearest_good_indices(single_phase_logph, dataset->coeffs_ph, _hmolar, _p, cached_single_phase_i,
911 cached_single_phase_j);
912 // Recalculate the phase
913 recalculate_singlephase_phase();
914 }
915 }
916 break;
917 }
918 case PT_INPUTS: {
919 _p = val1;
920 _T = val2;
921 if (!single_phase_logpT.native_inputs_are_in_range(_T, _p)) {
922 // Use the AbstractState instance
923 using_single_phase_table = false;
924 if (get_debug_level() > 5) {
925 std::cout << "inputs are not in range";
926 }
927 throw ValueError(format("inputs are not in range, p=%g Pa, T=%g K", _p, _T));
928 } else {
929 using_single_phase_table = true; // Use the table !
930 std::size_t iL = 0, iV = 0, iclosest = 0;
931 CoolPropDbl TL = 0, TV = 0;
932 SimpleState closest_state;
933 bool is_two_phase = false;
934 // Phase is imposed, use it
935 if (imposed_phase_index != iphase_not_imposed) {
936 is_two_phase = (imposed_phase_index == iphase_twophase);
937 } else {
938 if (is_mixture) {
939 is_two_phase = PhaseEnvelopeRoutines::is_inside(phase_envelope, iP, _p, iT, _T, iclosest, closest_state);
940 } else {
941 is_two_phase = pure_saturation.is_inside(iP, _p, iT, _T, iL, iV, TL, TV);
942 }
943 }
944 if (is_two_phase) {
945 using_single_phase_table = false;
946 throw ValueError(format("P,T with TTSE cannot be two-phase for now"));
947 } else {
948 selected_table = SELECTED_PT_TABLE;
949 // Find and cache the indices i, j
950 find_native_nearest_good_indices(single_phase_logpT, dataset->coeffs_pT, _T, _p, cached_single_phase_i, cached_single_phase_j);
951
952 if (imposed_phase_index != iphase_not_imposed) {
953 double rhoc = rhomolar_critical();
954 if (imposed_phase_index == iphase_liquid && cached_single_phase_i > 0) {
955 // We want a liquid solution, but we got a vapor solution
956 if (_p < this->AS->p_critical()) {
957 while (cached_single_phase_i > 0
958 && single_phase_logpT.rhomolar[cached_single_phase_i + 1][cached_single_phase_j] < rhoc) {
959 // Bump to lower temperature
960 cached_single_phase_i--;
961 }
962 double rho = evaluate_single_phase_pT(iDmolar, cached_single_phase_i, cached_single_phase_j);
963 if (rho < rhoc) {
964 // Didn't work
965 throw ValueError("Bump unsuccessful");
966 } else {
967 _rhomolar = rho;
968 }
969 }
970 } else if ((imposed_phase_index == iphase_gas || imposed_phase_index == iphase_supercritical_gas)
971 && cached_single_phase_i > 0) {
972
973 // We want a gas solution, but we got a liquid solution
974 if (_p < this->AS->p_critical()) {
975 while (cached_single_phase_i > 0
976 && single_phase_logpT.rhomolar[cached_single_phase_i][cached_single_phase_j + 1] > rhoc) {
977 // Bump to lower temperature
978 cached_single_phase_i++;
979 }
980 double rho = evaluate_single_phase_pT(iDmolar, cached_single_phase_i, cached_single_phase_j);
981 if (rho > rhoc) {
982 // Didn't work
983 throw ValueError("Bump unsuccessful");
984 } else {
985 _rhomolar = rho;
986 }
987 }
988 }
989 } else {
990
991 // If p < pc, you might be getting a liquid solution when you want a vapor solution or vice versa
992 // if you are very close to the saturation curve, so we figure out what the saturation temperature
993 // is for the given pressure
994 if (_p < this->AS->p_critical()) {
995 double Ts = pure_saturation.evaluate(iT, _p, _Q, iL, iV);
996 double TL = single_phase_logpT.T[cached_single_phase_i][cached_single_phase_j];
997 double TR = single_phase_logpT.T[cached_single_phase_i + 1][cached_single_phase_j];
998 if (TL < Ts && Ts < TR) {
999 if (_T < Ts) {
1000 if (cached_single_phase_i == 0) {
1001 throw ValueError(format("P, T are near saturation, but cannot move the cell to the left"));
1002 }
1003 // It's liquid, move the cell to the left
1004 cached_single_phase_i--;
1005 } else {
1006 if (cached_single_phase_i > single_phase_logpT.Nx - 2) {
1007 throw ValueError(format("P,T are near saturation, but cannot move the cell to the right"));
1008 }
1009 // It's vapor, move to the right
1010 cached_single_phase_i++;
1011 }
1012 }
1013 }
1014 }
1015 // Recalculate the phase
1016 recalculate_singlephase_phase();
1017 }
1018 }
1019 break;
1020 }
1021 case PUmolar_INPUTS:
1022 case PSmolar_INPUTS:
1023 case DmolarP_INPUTS: {
1024 CoolPropDbl otherval;
1025 parameters otherkey;
1026 switch (input_pair) {
1027 case PUmolar_INPUTS:
1028 _p = val1;
1029 _umolar = val2;
1030 otherval = val2;
1031 otherkey = iUmolar;
1032 break;
1033 case PSmolar_INPUTS:
1034 _p = val1;
1035 _smolar = val2;
1036 otherval = val2;
1037 otherkey = iSmolar;
1038 break;
1039 case DmolarP_INPUTS:
1040 _rhomolar = val1;
1041 _p = val2;
1042 otherval = val1;
1043 otherkey = iDmolar;
1044 break;
1045 default:
1046 throw ValueError("Bad (impossible) pair");
1047 }
1048
1049 using_single_phase_table = true; // Use the table (or first guess is that it is single-phase)!
1050 std::size_t iL = std::numeric_limits<std::size_t>::max(), iV = std::numeric_limits<std::size_t>::max();
1051 CoolPropDbl zL = 0, zV = 0;
1052 std::size_t iclosest = 0;
1053 SimpleState closest_state;
1054 bool is_two_phase = false;
1055 // If phase is imposed, use it, but only if it's single phase:
1056 // - Imposed two phase still needs to determine saturation limits by calling is_inside().
1057 // - There's no speed increase to be gained by imposing two phase.
1058 if ((imposed_phase_index == iphase_not_imposed) || (imposed_phase_index == iphase_twophase)) {
1059 if (is_mixture) {
1060 is_two_phase = PhaseEnvelopeRoutines::is_inside(phase_envelope, iP, _p, otherkey, otherval, iclosest, closest_state);
1061 if (is_two_phase) {
1062 std::vector<std::pair<std::size_t, std::size_t>> intersect =
1063 PhaseEnvelopeRoutines::find_intersections(phase_envelope, iP, _p);
1064 if (intersect.size() < 2) {
1065 throw ValueError(format("p [%g Pa] is not within phase envelope", _p));
1066 }
1067 iV = intersect[0].first;
1068 iL = intersect[1].first;
1069 zL = PhaseEnvelopeRoutines::evaluate(phase_envelope, otherkey, iP, _p, iL);
1070 zV = PhaseEnvelopeRoutines::evaluate(phase_envelope, otherkey, iP, _p, iV);
1071 }
1072 } else {
1073 is_two_phase = pure_saturation.is_inside(iP, _p, otherkey, otherval, iL, iV, zL, zV);
1074 }
1075 }
1076 // Phase determined or imposed, now interpolate results
1077 if (is_two_phase) {
1078 using_single_phase_table = false;
1079 if (otherkey == iDmolar) {
1080 _Q = (1 / otherval - 1 / zL) / (1 / zV - 1 / zL);
1081 } else {
1082 _Q = (otherval - zL) / (zV - zL);
1083 }
1084 if (!is_in_closed_range(0.0, 1.0, static_cast<double>(_Q))) {
1085
1086 throw ValueError(format("vapor quality is not in (0,1) for %s: %g p: %g", get_parameter_information(otherkey, "short").c_str(),
1087 otherval, static_cast<double>(_p)));
1088 } else if (!is_mixture) {
1089 cached_saturation_iL = iL;
1090 cached_saturation_iV = iV;
1091 }
1092 _phase = iphase_twophase;
1093 } else {
1094 selected_table = SELECTED_PH_TABLE;
1095 // Find and cache the indices i, j
1096 find_nearest_neighbor(single_phase_logph, dataset->coeffs_ph, iP, _p, otherkey, otherval, cached_single_phase_i,
1097 cached_single_phase_j);
1098 // Now find hmolar given P, X for X in Smolar, Umolar, Dmolar
1099 invert_single_phase_x(single_phase_logph, dataset->coeffs_ph, otherkey, otherval, _p, cached_single_phase_i, cached_single_phase_j);
1100 // Recalculate the phase
1101 recalculate_singlephase_phase();
1102 }
1103 break;
1104 }
1105 case SmolarT_INPUTS:
1106 case DmolarT_INPUTS: {
1107 CoolPropDbl otherval;
1108 parameters otherkey;
1109 switch (input_pair) {
1110 case SmolarT_INPUTS:
1111 _smolar = val1;
1112 _T = val2;
1113 otherval = val1;
1114 otherkey = iSmolar;
1115 break;
1116 case DmolarT_INPUTS:
1117 _rhomolar = val1;
1118 _T = val2;
1119 otherval = val1;
1120 otherkey = iDmolar;
1121 break;
1122 default:
1123 throw ValueError("Bad (impossible) pair");
1124 }
1125
1126 using_single_phase_table = true; // Use the table (or first guess is that it is single-phase)!
1127 std::size_t iL = std::numeric_limits<std::size_t>::max(), iV = std::numeric_limits<std::size_t>::max();
1128 CoolPropDbl zL = 0, zV = 0;
1129 std::size_t iclosest = 0;
1130 SimpleState closest_state;
1131 bool is_two_phase = false;
1132 // If phase is imposed, use it, but only if it's single phase:
1133 // - Imposed two phase still needs to determine saturation limits by calling is_inside().
1134 // - There's no speed increase to be gained by imposing two phase.
1135 if ((imposed_phase_index == iphase_not_imposed) || (imposed_phase_index == iphase_twophase)) {
1136 if (is_mixture) {
1137 is_two_phase = PhaseEnvelopeRoutines::is_inside(phase_envelope, iT, _T, otherkey, otherval, iclosest, closest_state);
1138 if (is_two_phase) {
1139 std::vector<std::pair<std::size_t, std::size_t>> intersect =
1140 PhaseEnvelopeRoutines::find_intersections(phase_envelope, iT, _T);
1141 if (intersect.size() < 2) {
1142 throw ValueError(format("T [%g K] is not within phase envelope", _T));
1143 }
1144 iV = intersect[0].first;
1145 iL = intersect[1].first;
1146 zL = PhaseEnvelopeRoutines::evaluate(phase_envelope, otherkey, iT, _T, iL);
1147 zV = PhaseEnvelopeRoutines::evaluate(phase_envelope, otherkey, iT, _T, iV);
1148 }
1149 } else {
1150 is_two_phase = pure_saturation.is_inside(iT, _T, otherkey, otherval, iL, iV, zL, zV);
1151 }
1152 }
1153 if (is_two_phase) {
1154 using_single_phase_table = false;
1155 if (otherkey == iDmolar) {
1156 _Q = (1 / otherval - 1 / zL) / (1 / zV - 1 / zL);
1157 } else {
1158 _Q = (otherval - zL) / (zV - zL);
1159 }
1160 if (!is_in_closed_range(0.0, 1.0, static_cast<double>(_Q))) {
1161 throw ValueError(format("vapor quality is not in (0,1) for %s: %g T: %g", get_parameter_information(otherkey, "short").c_str(),
1162 otherval, static_cast<double>(_T)));
1163 } else if (!is_mixture) {
1164 cached_saturation_iL = iL;
1165 cached_saturation_iV = iV;
1166 _p = pure_saturation.evaluate(iP, _T, _Q, iL, iV);
1167 } else {
1168 // Mixture
1169 std::vector<std::pair<std::size_t, std::size_t>> intersect = PhaseEnvelopeRoutines::find_intersections(phase_envelope, iT, _T);
1170 if (intersect.size() < 2) {
1171 throw ValueError(format("T [%g K] is not within phase envelope", _T));
1172 }
1173 iV = intersect[0].first;
1174 iL = intersect[1].first;
1175 CoolPropDbl pL = PhaseEnvelopeRoutines::evaluate(phase_envelope, iP, iT, _T, iL);
1176 CoolPropDbl pV = PhaseEnvelopeRoutines::evaluate(phase_envelope, iP, iT, _T, iV);
1177 _p = _Q * pV + (1 - _Q) * pL;
1178 }
1179 } else {
1180 selected_table = SELECTED_PT_TABLE;
1181 // Find and cache the indices i, j
1182 find_nearest_neighbor(single_phase_logpT, dataset->coeffs_pT, iT, _T, otherkey, otherval, cached_single_phase_i,
1183 cached_single_phase_j);
1184 // Now find the y variable (Dmolar or Smolar in this case)
1185 invert_single_phase_y(single_phase_logpT, dataset->coeffs_pT, otherkey, otherval, _T, cached_single_phase_i, cached_single_phase_j);
1186 // Recalculate the phase
1187 recalculate_singlephase_phase();
1188 }
1189 break;
1190 }
1191 case PQ_INPUTS: {
1192 std::size_t iL = 0, iV = 0;
1193 _p = val1;
1194 _Q = val2;
1195 using_single_phase_table = false;
1196 if (!is_in_closed_range(0.0, 1.0, static_cast<double>(_Q))) {
1197 throw ValueError(format("vapor quality [%g] is not in (0,1)", static_cast<double>(val2)));
1198 } else {
1199 CoolPropDbl TL = _HUGE, TV = _HUGE;
1200 if (is_mixture) {
1201 std::vector<std::pair<std::size_t, std::size_t>> intersect = PhaseEnvelopeRoutines::find_intersections(phase_envelope, iP, _p);
1202 if (intersect.size() < 2) {
1203 throw ValueError(format("p [%g Pa] is not within phase envelope", _p));
1204 }
1205 iV = intersect[0].first;
1206 iL = intersect[1].first;
1207 TL = PhaseEnvelopeRoutines::evaluate(phase_envelope, iT, iP, _p, iL);
1208 TV = PhaseEnvelopeRoutines::evaluate(phase_envelope, iT, iP, _p, iV);
1209 } else {
1210 bool it_is_inside = pure_saturation.is_inside(iP, _p, iQ, _Q, iL, iV, TL, TV);
1211 if (!it_is_inside) {
1212 throw ValueError("Not possible to determine whether pressure is inside or not");
1213 }
1214 }
1215 _T = _Q * TV + (1 - _Q) * TL;
1216 cached_saturation_iL = iL;
1217 cached_saturation_iV = iV;
1218 _phase = iphase_twophase;
1219 }
1220 break;
1221 }
1222 case QT_INPUTS: {
1223 std::size_t iL = 0, iV = 0;
1224 _Q = val1;
1225 _T = val2;
1226
1227 using_single_phase_table = false;
1228 if (!is_in_closed_range(0.0, 1.0, static_cast<double>(_Q))) {
1229 throw ValueError(format("vapor quality [%g] is not in (0,1)", static_cast<double>(val1)));
1230 } else {
1231 CoolPropDbl pL = _HUGE, pV = _HUGE;
1232 if (is_mixture) {
1233 std::vector<std::pair<std::size_t, std::size_t>> intersect = PhaseEnvelopeRoutines::find_intersections(phase_envelope, iT, _T);
1234 if (intersect.size() < 2) {
1235 throw ValueError(format("T [%g K] is not within phase envelope", _T));
1236 }
1237 iV = intersect[0].first;
1238 iL = intersect[1].first;
1239 pL = PhaseEnvelopeRoutines::evaluate(phase_envelope, iP, iT, _T, iL);
1240 pV = PhaseEnvelopeRoutines::evaluate(phase_envelope, iP, iT, _T, iV);
1241 } else {
1242 pure_saturation.is_inside(iT, _T, iQ, _Q, iL, iV, pL, pV);
1243 }
1244 _p = _Q * pV + (1 - _Q) * pL;
1245 cached_saturation_iL = iL;
1246 cached_saturation_iV = iV;
1247 _phase = iphase_twophase;
1248 }
1249 break;
1250 }
1251 default:
1252 throw ValueError("Sorry, but this set of inputs is not supported for Tabular backend");
1253 }
1254}
1255
1256void CoolProp::TabularDataSet::write_tables(const std::string& path_to_tables) {
1257 make_dirs(path_to_tables);
1258 write_table(single_phase_logph, path_to_tables, "single_phase_logph");
1259 write_table(single_phase_logpT, path_to_tables, "single_phase_logpT");
1260 write_table(pure_saturation, path_to_tables, "pure_saturation");
1261 write_table(phase_envelope, path_to_tables, "phase_envelope");
1262}
1263
1264void CoolProp::TabularDataSet::load_tables(const std::string& path_to_tables, shared_ptr<CoolProp::AbstractState>& AS) {
1265 single_phase_logph.AS = AS;
1266 single_phase_logpT.AS = AS;
1267 pure_saturation.AS = AS;
1268 single_phase_logph.set_limits();
1269 single_phase_logpT.set_limits();
1270 load_table(single_phase_logph, path_to_tables, "single_phase_logph.bin.z");
1271 load_table(single_phase_logpT, path_to_tables, "single_phase_logpT.bin.z");
1272 load_table(pure_saturation, path_to_tables, "pure_saturation.bin.z");
1273 load_table(phase_envelope, path_to_tables, "phase_envelope.bin.z");
1274 tables_loaded = true;
1275 if (get_debug_level() > 0) {
1276 std::cout << "Tables loaded" << std::endl;
1277 }
1278};
1279
1280void CoolProp::TabularDataSet::build_tables(shared_ptr<CoolProp::AbstractState>& AS) {
1281 // Pure or pseudo-pure fluid
1282 if (AS->get_mole_fractions().size() == 1) {
1283 pure_saturation.build(AS);
1284 } else {
1285 // Call function to actually construct the phase envelope
1286 AS->build_phase_envelope("");
1287 // Copy constructed phase envelope into this class
1288 PhaseEnvelopeData PED = AS->get_phase_envelope_data();
1289 // Convert into packable form
1290 phase_envelope.copy_from_nonpackable(PED);
1291 // Resize so that it will load properly
1292 pure_saturation.resize(pure_saturation.N);
1293 }
1294 single_phase_logph.build(AS);
1295 single_phase_logpT.build(AS);
1296 tables_loaded = true;
1297}
1298
1300CoolProp::TabularDataSet* CoolProp::TabularDataLibrary::get_set_of_tables(shared_ptr<AbstractState>& AS, bool& loaded) {
1301 const std::string path = path_to_tables(AS);
1302 // Try to find tabular set if it is already loaded
1303 std::map<std::string, TabularDataSet>::iterator it = data.find(path);
1304 // It is already in the map, return it
1305 if (it != data.end()) {
1306 loaded = it->second.tables_loaded;
1307 return &(it->second);
1308 }
1309 // It is not in the map, build it
1310 else {
1311 TabularDataSet set;
1312 data.insert(std::pair<std::string, TabularDataSet>(path, set));
1313 TabularDataSet& dataset = data[path];
1314 try {
1315 if (!dataset.tables_loaded) {
1316 dataset.load_tables(path, AS);
1317 }
1318 loaded = true;
1319 } catch (std::exception&) {
1320 loaded = false;
1321 }
1322 return &(dataset);
1323 }
1324}
1325
1326void CoolProp::TabularDataSet::build_coeffs(SinglePhaseGriddedTableData& table, std::vector<std::vector<CellCoeffs>>& coeffs) {
1327 if (!coeffs.empty()) {
1328 return;
1329 }
1330 const bool debug = get_debug_level() > 5 || false;
1331 const int param_count = 6;
1332 parameters param_list[param_count] = {iDmolar, iT, iSmolar, iHmolar, iP, iUmolar};
1333 std::vector<std::vector<double>>*f = NULL, *fx = NULL, *fy = NULL, *fxy = NULL;
1334
1335 clock_t t1 = clock();
1336
1337 // Resize the coefficient structures
1338 coeffs.resize(table.Nx - 1, std::vector<CellCoeffs>(table.Ny - 1));
1339
1340 int valid_cell_count = 0;
1341 for (std::size_t k = 0; k < param_count; ++k) {
1342 parameters param = param_list[k];
1343 if (param == table.xkey || param == table.ykey) {
1344 continue;
1345 } // Skip tables that match either of the input variables
1346
1347 switch (param) {
1348 case iT:
1349 f = &(table.T);
1350 fx = &(table.dTdx);
1351 fy = &(table.dTdy);
1352 fxy = &(table.d2Tdxdy);
1353 break;
1354 case iP:
1355 f = &(table.p);
1356 fx = &(table.dpdx);
1357 fy = &(table.dpdy);
1358 fxy = &(table.d2pdxdy);
1359 break;
1360 case iDmolar:
1361 f = &(table.rhomolar);
1362 fx = &(table.drhomolardx);
1363 fy = &(table.drhomolardy);
1364 fxy = &(table.d2rhomolardxdy);
1365 break;
1366 case iSmolar:
1367 f = &(table.smolar);
1368 fx = &(table.dsmolardx);
1369 fy = &(table.dsmolardy);
1370 fxy = &(table.d2smolardxdy);
1371 break;
1372 case iHmolar:
1373 f = &(table.hmolar);
1374 fx = &(table.dhmolardx);
1375 fy = &(table.dhmolardy);
1376 fxy = &(table.d2hmolardxdy);
1377 break;
1378 case iUmolar:
1379 f = &(table.umolar);
1380 fx = &(table.dumolardx);
1381 fy = &(table.dumolardy);
1382 fxy = &(table.d2umolardxdy);
1383 break;
1384 default:
1385 throw ValueError("Invalid variable type to build_coeffs");
1386 }
1387 for (std::size_t i = 0; i < table.Nx - 1; ++i) // -1 since we have one fewer cells than nodes
1388 {
1389 for (std::size_t j = 0; j < table.Ny - 1; ++j) // -1 since we have one fewer cells than nodes
1390 {
1391 if (ValidNumber((*f)[i][j]) && ValidNumber((*f)[i + 1][j]) && ValidNumber((*f)[i][j + 1]) && ValidNumber((*f)[i + 1][j + 1])) {
1392
1393 // This will hold the scaled f values for the cell
1394 Eigen::Matrix<double, 16, 1> F;
1395 // The output values (do not require scaling
1396 F(0) = (*f)[i][j];
1397 F(1) = (*f)[i + 1][j];
1398 F(2) = (*f)[i][j + 1];
1399 F(3) = (*f)[i + 1][j + 1];
1400 // Scaling parameter
1401 // d(f)/dxhat = df/dx * dx/dxhat, where xhat = (x-x_i)/(x_{i+1}-x_i)
1402 coeffs[i][j].dx_dxhat = table.xvec[i + 1] - table.xvec[i];
1403 double dx_dxhat = coeffs[i][j].dx_dxhat;
1404 F(4) = (*fx)[i][j] * dx_dxhat;
1405 F(5) = (*fx)[i + 1][j] * dx_dxhat;
1406 F(6) = (*fx)[i][j + 1] * dx_dxhat;
1407 F(7) = (*fx)[i + 1][j + 1] * dx_dxhat;
1408 // Scaling parameter
1409 // d(f)/dyhat = df/dy * dy/dyhat, where yhat = (y-y_j)/(y_{j+1}-y_j)
1410 coeffs[i][j].dy_dyhat = table.yvec[j + 1] - table.yvec[j];
1411 double dy_dyhat = coeffs[i][j].dy_dyhat;
1412 F(8) = (*fy)[i][j] * dy_dyhat;
1413 F(9) = (*fy)[i + 1][j] * dy_dyhat;
1414 F(10) = (*fy)[i][j + 1] * dy_dyhat;
1415 F(11) = (*fy)[i + 1][j + 1] * dy_dyhat;
1416 // Cross derivatives are doubly scaled following the examples above
1417 F(12) = (*fxy)[i][j] * dy_dyhat * dx_dxhat;
1418 F(13) = (*fxy)[i + 1][j] * dy_dyhat * dx_dxhat;
1419 F(14) = (*fxy)[i][j + 1] * dy_dyhat * dx_dxhat;
1420 F(15) = (*fxy)[i + 1][j + 1] * dy_dyhat * dx_dxhat;
1421 // Calculate the alpha coefficients
1422 Eigen::MatrixXd alpha = Ainv.transpose() * F; // 16x1; Watch out for the transpose!
1423 std::vector<double> valpha = eigen_to_vec1D(alpha);
1424 coeffs[i][j].set(param, valpha);
1425 coeffs[i][j].set_valid();
1426 valid_cell_count++;
1427 } else {
1428 coeffs[i][j].set_invalid();
1429 }
1430 }
1431 }
1432 double elapsed = (clock() - t1) / ((double)CLOCKS_PER_SEC);
1433 if (debug) {
1434 std::cout << format("Calculated bicubic coefficients for %d good cells in %g sec.\n", valid_cell_count, elapsed);
1435 }
1436 std::size_t remap_count = 0;
1437 // Now find invalid cells and give them pointers to a neighboring cell that works
1438 for (std::size_t i = 0; i < table.Nx - 1; ++i) // -1 since we have one fewer cells than nodes
1439 {
1440 for (std::size_t j = 0; j < table.Ny - 1; ++j) // -1 since we have one fewer cells than nodes
1441 {
1442 // Not a valid cell
1443 if (!coeffs[i][j].valid()) {
1444 // Offsets that we are going to try in order (left, right, top, bottom, diagonals)
1445 int xoffsets[] = {-1, 1, 0, 0, -1, 1, 1, -1};
1446 int yoffsets[] = {0, 0, 1, -1, -1, -1, 1, 1};
1447 // Length of offset
1448 std::size_t N = sizeof(xoffsets) / sizeof(xoffsets[0]);
1449 for (std::size_t k = 0; k < N; ++k) {
1450 std::size_t iplus = i + xoffsets[k];
1451 std::size_t jplus = j + yoffsets[k];
1452 if (0 < iplus && iplus < table.Nx - 1 && 0 < jplus && jplus < table.Ny - 1 && coeffs[iplus][jplus].valid()) {
1453 coeffs[i][j].set_alternate(iplus, jplus);
1454 remap_count++;
1455 if (debug) {
1456 std::cout << format("Mapping %d,%d to %d,%d\n", i, j, iplus, jplus);
1457 }
1458 break;
1459 }
1460 }
1461 }
1462 }
1463 }
1464 if (debug) {
1465 std::cout << format("Remapped %d cells\n", remap_count);
1466 }
1467 }
1468}
1469
1470# if defined(ENABLE_CATCH)
1471# include <catch2/catch_all.hpp>
1472
1473// Defined global so we only load once
1474static shared_ptr<CoolProp::AbstractState> ASHEOS, ASTTSE, ASBICUBIC;
1475
1476/* Use a fixture so that loading of the tables from memory only happens once in the initializer */
1477class TabularFixture
1478{
1479 public:
1480 TabularFixture() {}
1481 void setup() {
1482 if (ASHEOS.get() == NULL) {
1483 ASHEOS.reset(CoolProp::AbstractState::factory("HEOS", "Water"));
1484 }
1485 if (ASTTSE.get() == NULL) {
1486 ASTTSE.reset(CoolProp::AbstractState::factory("TTSE&HEOS", "Water"));
1487 }
1488 if (ASBICUBIC.get() == NULL) {
1489 ASBICUBIC.reset(CoolProp::AbstractState::factory("BICUBIC&HEOS", "Water"));
1490 }
1491 }
1492};
1493TEST_CASE_METHOD(TabularFixture, "Tests for tabular backends with water", "[Tabular]") {
1494 SECTION("first_saturation_deriv invalid quality") {
1495 setup();
1496 ASBICUBIC->update(CoolProp::PQ_INPUTS, 101325, 0.5);
1497 CHECK_THROWS(ASBICUBIC->first_saturation_deriv(CoolProp::iP, CoolProp::iT));
1498 ASTTSE->update(CoolProp::PQ_INPUTS, 101325, 0.5);
1499 CHECK_THROWS(ASTTSE->first_saturation_deriv(CoolProp::iP, CoolProp::iT));
1500 }
1501
1502 SECTION("first_saturation_deriv dp/dT") {
1503 setup();
1504 ASHEOS->update(CoolProp::PQ_INPUTS, 101325, 1);
1505 CoolPropDbl expected = ASHEOS->first_saturation_deriv(CoolProp::iP, CoolProp::iT);
1506 ASTTSE->update(CoolProp::PQ_INPUTS, 101325, 1);
1507 CoolPropDbl actual_TTSE = ASTTSE->first_saturation_deriv(CoolProp::iP, CoolProp::iT);
1508 ASTTSE->update(CoolProp::PQ_INPUTS, 101325, 1);
1509 CoolPropDbl actual_BICUBIC = ASTTSE->first_saturation_deriv(CoolProp::iP, CoolProp::iT);
1510 CAPTURE(expected);
1511 CAPTURE(actual_TTSE);
1512 CAPTURE(actual_BICUBIC);
1513 CHECK(std::abs((expected - actual_TTSE) / expected) < 1e-6);
1514 CHECK(std::abs((expected - actual_BICUBIC) / expected) < 1e-6);
1515 }
1516 SECTION("first_saturation_deriv dDmolar/dP") {
1517 setup();
1518 ASHEOS->update(CoolProp::PQ_INPUTS, 101325, 1);
1519 CoolPropDbl expected = ASHEOS->first_saturation_deriv(CoolProp::iDmolar, CoolProp::iP);
1520 ASTTSE->update(CoolProp::PQ_INPUTS, 101325, 1);
1521 CoolPropDbl actual_TTSE = ASTTSE->first_saturation_deriv(CoolProp::iDmolar, CoolProp::iP);
1522 ASBICUBIC->update(CoolProp::PQ_INPUTS, 101325, 1);
1523 CoolPropDbl actual_BICUBIC = ASBICUBIC->first_saturation_deriv(CoolProp::iDmolar, CoolProp::iP);
1524 CAPTURE(expected);
1525 CAPTURE(actual_TTSE);
1526 CAPTURE(actual_BICUBIC);
1527 CHECK(std::abs((expected - actual_TTSE) / expected) < 1e-6);
1528 CHECK(std::abs((expected - actual_BICUBIC) / expected) < 1e-6);
1529 }
1530 SECTION("first_saturation_deriv dDmass/dP") {
1531 setup();
1532 ASHEOS->update(CoolProp::PQ_INPUTS, 101325, 1);
1533 CoolPropDbl expected = ASHEOS->first_saturation_deriv(CoolProp::iDmass, CoolProp::iP);
1534 ASTTSE->update(CoolProp::PQ_INPUTS, 101325, 1);
1535 CoolPropDbl actual_TTSE = ASTTSE->first_saturation_deriv(CoolProp::iDmass, CoolProp::iP);
1536 ASBICUBIC->update(CoolProp::PQ_INPUTS, 101325, 1);
1537 CoolPropDbl actual_BICUBIC = ASBICUBIC->first_saturation_deriv(CoolProp::iDmass, CoolProp::iP);
1538 CAPTURE(expected);
1539 CAPTURE(actual_TTSE);
1540 CAPTURE(actual_BICUBIC);
1541 CHECK(std::abs((expected - actual_TTSE) / expected) < 1e-6);
1542 CHECK(std::abs((expected - actual_BICUBIC) / expected) < 1e-6);
1543 }
1544 SECTION("first_saturation_deriv dHmass/dP") {
1545 setup();
1546 ASHEOS->update(CoolProp::PQ_INPUTS, 101325, 1);
1547 CoolPropDbl expected = ASHEOS->first_saturation_deriv(CoolProp::iHmass, CoolProp::iP);
1548 ASTTSE->update(CoolProp::PQ_INPUTS, 101325, 1);
1549 CoolPropDbl actual_TTSE = ASTTSE->first_saturation_deriv(CoolProp::iHmass, CoolProp::iP);
1550 ASBICUBIC->update(CoolProp::PQ_INPUTS, 101325, 1);
1551 CoolPropDbl actual_BICUBIC = ASBICUBIC->first_saturation_deriv(CoolProp::iHmass, CoolProp::iP);
1552 CAPTURE(expected);
1553 CAPTURE(actual_TTSE);
1554 CAPTURE(actual_BICUBIC);
1555 CHECK(std::abs((expected - actual_TTSE) / expected) < 1e-6);
1556 CHECK(std::abs((expected - actual_BICUBIC) / expected) < 1e-6);
1557 }
1558 SECTION("first_saturation_deriv dHmass/dP w/ QT as inputs") {
1559 setup();
1560 ASHEOS->update(CoolProp::QT_INPUTS, 1, 370);
1561 CoolPropDbl expected = ASHEOS->first_saturation_deriv(CoolProp::iHmass, CoolProp::iP);
1562 ASTTSE->update(CoolProp::QT_INPUTS, 1, 370);
1563 CoolPropDbl actual_TTSE = ASTTSE->first_saturation_deriv(CoolProp::iHmass, CoolProp::iP);
1564 ASBICUBIC->update(CoolProp::QT_INPUTS, 1, 370);
1565 CoolPropDbl actual_BICUBIC = ASBICUBIC->first_saturation_deriv(CoolProp::iHmass, CoolProp::iP);
1566 CAPTURE(expected);
1567 CAPTURE(actual_TTSE);
1568 CAPTURE(actual_BICUBIC);
1569 CHECK(std::abs((expected - actual_TTSE) / expected) < 1e-6);
1570 CHECK(std::abs((expected - actual_BICUBIC) / expected) < 1e-6);
1571 }
1572 SECTION("first_two_phase_deriv dDmolar/dP|Hmolar") {
1573 setup();
1574 ASHEOS->update(CoolProp::PQ_INPUTS, 101325, 0.1);
1575 CoolPropDbl expected = ASHEOS->first_two_phase_deriv(CoolProp::iDmolar, CoolProp::iP, CoolProp::iHmolar);
1576 ASTTSE->update(CoolProp::PQ_INPUTS, 101325, 0.1);
1577 CoolPropDbl actual_TTSE = ASTTSE->first_two_phase_deriv(CoolProp::iDmolar, CoolProp::iP, CoolProp::iHmolar);
1578 ASBICUBIC->update(CoolProp::PQ_INPUTS, 101325, 0.1);
1579 CoolPropDbl actual_BICUBIC = ASBICUBIC->first_two_phase_deriv(CoolProp::iDmolar, CoolProp::iP, CoolProp::iHmolar);
1580 CAPTURE(expected);
1581 CAPTURE(actual_TTSE);
1582 CAPTURE(actual_BICUBIC);
1583 CHECK(std::abs((expected - actual_TTSE) / expected) < 1e-6);
1584 CHECK(std::abs((expected - actual_BICUBIC) / expected) < 1e-6);
1585 }
1586 SECTION("first_two_phase_deriv dDmass/dP|Hmass") {
1587 setup();
1588 ASHEOS->update(CoolProp::PQ_INPUTS, 101325, 0.1);
1589 CoolPropDbl expected = ASHEOS->first_two_phase_deriv(CoolProp::iDmass, CoolProp::iP, CoolProp::iHmass);
1590 ASTTSE->update(CoolProp::PQ_INPUTS, 101325, 0.1);
1591 CoolPropDbl actual_TTSE = ASTTSE->first_two_phase_deriv(CoolProp::iDmass, CoolProp::iP, CoolProp::iHmass);
1592 ASBICUBIC->update(CoolProp::PQ_INPUTS, 101325, 0.1);
1593 CoolPropDbl actual_BICUBIC = ASBICUBIC->first_two_phase_deriv(CoolProp::iDmass, CoolProp::iP, CoolProp::iHmass);
1594 CAPTURE(expected);
1595 CAPTURE(actual_TTSE);
1596 CAPTURE(actual_BICUBIC);
1597 CHECK(std::abs((expected - actual_TTSE) / expected) < 1e-6);
1598 CHECK(std::abs((expected - actual_BICUBIC) / expected) < 1e-6);
1599 }
1600 SECTION("first_two_phase_deriv dDmass/dHmass|P") {
1601 setup();
1602 ASHEOS->update(CoolProp::PQ_INPUTS, 101325, 0.1);
1603 CoolPropDbl expected = ASHEOS->first_two_phase_deriv(CoolProp::iDmass, CoolProp::iHmass, CoolProp::iP);
1604 ASTTSE->update(CoolProp::PQ_INPUTS, 101325, 0.1);
1605 CoolPropDbl actual_TTSE = ASTTSE->first_two_phase_deriv(CoolProp::iDmass, CoolProp::iHmass, CoolProp::iP);
1606 ASBICUBIC->update(CoolProp::PQ_INPUTS, 101325, 0.1);
1607 CoolPropDbl actual_BICUBIC = ASBICUBIC->first_two_phase_deriv(CoolProp::iDmass, CoolProp::iHmass, CoolProp::iP);
1608 CAPTURE(expected);
1609 CAPTURE(actual_TTSE);
1610 CAPTURE(actual_BICUBIC);
1611 CHECK(std::abs((expected - actual_TTSE) / expected) < 1e-6);
1612 CHECK(std::abs((expected - actual_BICUBIC) / expected) < 1e-6);
1613 }
1614 SECTION("first_partial_deriv dHmass/dT|P") {
1615 setup();
1616 ASHEOS->update(CoolProp::PT_INPUTS, 101325, 300);
1617 double expected = ASHEOS->cpmass();
1618 ASTTSE->update(CoolProp::PT_INPUTS, 101325, 300);
1619 double dhdT_TTSE = ASTTSE->first_partial_deriv(CoolProp::iHmass, CoolProp::iT, CoolProp::iP);
1620 ASBICUBIC->update(CoolProp::PT_INPUTS, 101325, 300);
1621 double dhdT_BICUBIC = ASBICUBIC->first_partial_deriv(CoolProp::iHmass, CoolProp::iT, CoolProp::iP);
1622 CAPTURE(expected);
1623 CAPTURE(dhdT_TTSE);
1624 CAPTURE(dhdT_BICUBIC);
1625 CHECK(std::abs((expected - dhdT_TTSE) / expected) < 1e-4);
1626 CHECK(std::abs((expected - dhdT_BICUBIC) / expected) < 1e-4);
1627 }
1628 SECTION("first_partial_deriv dHmolar/dT|P") {
1629 setup();
1630 ASHEOS->update(CoolProp::PT_INPUTS, 101325, 300);
1631 double expected = ASHEOS->cpmolar();
1632 ASTTSE->update(CoolProp::PT_INPUTS, 101325, 300);
1633 double dhdT_TTSE = ASTTSE->first_partial_deriv(CoolProp::iHmolar, CoolProp::iT, CoolProp::iP);
1634 ASBICUBIC->update(CoolProp::PT_INPUTS, 101325, 300);
1635 double dhdT_BICUBIC = ASBICUBIC->first_partial_deriv(CoolProp::iHmolar, CoolProp::iT, CoolProp::iP);
1636 CAPTURE(expected);
1637 CAPTURE(dhdT_TTSE);
1638 CAPTURE(dhdT_BICUBIC);
1639 CHECK(std::abs((expected - dhdT_TTSE) / expected) < 1e-4);
1640 CHECK(std::abs((expected - dhdT_BICUBIC) / expected) < 1e-4);
1641 }
1642 SECTION("check isentropic process") {
1643 setup();
1644 double T0 = 300;
1645 double p0 = 1e5;
1646 double p1 = 1e6;
1647 ASHEOS->update(CoolProp::PT_INPUTS, p0, 300);
1648 double s0 = ASHEOS->smolar();
1649 ASHEOS->update(CoolProp::PSmolar_INPUTS, p1, s0);
1650 double expected = ASHEOS->T();
1651 ASTTSE->update(CoolProp::PSmolar_INPUTS, p1, s0);
1652 double actual_TTSE = ASTTSE->T();
1653 ASBICUBIC->update(CoolProp::PSmolar_INPUTS, p1, s0);
1654 double actual_BICUBIC = ASBICUBIC->T();
1655 CAPTURE(expected);
1656 CAPTURE(actual_TTSE);
1657 CAPTURE(actual_BICUBIC);
1658 CHECK(std::abs((expected - actual_TTSE) / expected) < 1e-2);
1659 CHECK(std::abs((expected - actual_BICUBIC) / expected) < 1e-2);
1660 }
1661 SECTION("check D=1 mol/m3, T=500 K inputs") {
1662 setup();
1663 double d = 1;
1664 CAPTURE(d);
1665 ASHEOS->update(CoolProp::DmolarT_INPUTS, d, 500);
1666 double expected = ASHEOS->p();
1667 ASTTSE->update(CoolProp::DmolarT_INPUTS, d, 500);
1668 double actual_TTSE = ASTTSE->p();
1669 ASBICUBIC->update(CoolProp::DmolarT_INPUTS, d, 500);
1670 double actual_BICUBIC = ASBICUBIC->p();
1671 CAPTURE(expected);
1672 CAPTURE(actual_TTSE);
1673 CAPTURE(actual_BICUBIC);
1674 CHECK(std::abs((expected - actual_TTSE) / expected) < 1e-3);
1675 CHECK(std::abs((expected - actual_BICUBIC) / expected) < 1e-3);
1676 }
1677}
1678# endif // ENABLE_CATCH
1679
1680#endif // !defined(NO_TABULAR_BACKENDS)
Generated on Mon Jun 2 2025 05:21:52 for CoolProp by doxygen 1.9.4