CoolProp 8.0.0
An open-source fluid property and humid air property database
IF97Backend.h
Go to the documentation of this file.
1
2#ifndef IF97BACKEND_H_
3#define IF97BACKEND_H_
4
5#include "CoolProp/DataStructures.h"
6#include <IF97.h>
7#include "CoolProp/AbstractState.h"
8#include "CoolProp/Exceptions.h"
9#include <vector>
10
11namespace CoolProp {
12
13class IF97Backend : public AbstractState
14{
15
16 protected:
19 CachedElement _hmass, _rhomass, _smass,
20 _reverse; // Also need a way to flag using the IF97 reverse calcs for h(p,s) and s(p,h)
22
23 public:
25 std::string backend_name() override {
26 return get_backend_string(IF97_BACKEND);
27 }
28
34 std::vector<std::string> calc_fluid_names() override {
35 return {"Water"};
36 }
37
38 // REQUIRED BUT NOT USED IN IF97 FUNCTIONS
39 bool using_mole_fractions() override {
40 return false;
41 };
42 bool using_mass_fractions() override {
43 return true;
44 }; // But actually it doesn't matter since it is only a pure fluid
45 bool using_volu_fractions() override {
46 return false;
47 };
48 void set_mole_fractions(const std::vector<CoolPropDbl>& mole_fractions) override {
49 throw NotImplementedError("Mole composition has not been implemented.");
50 };
51 void set_mass_fractions(const std::vector<CoolPropDbl>& mass_fractions) override {}; // Not implemented, but don't throw any errors
52 void set_volu_fractions(const std::vector<CoolPropDbl>& volu_fractions) override {
53 throw NotImplementedError("Volume composition has not been implemented.");
54 };
55 const std::vector<CoolPropDbl>& get_mole_fractions() override {
56 throw NotImplementedError("get_mole_fractions composition has not been implemented.");
57 };
58
60 bool clear() override {
61 // Reset the full CachedElement registry so every cached value
62 // (speed_sound, cpmass, cvmass, umass, ...) gets invalidated
63 // alongside the IF97-specific fields the override knows about.
64 // Without `cache.clear()`, a backend instance reused across
65 // multiple update() calls would return stale property values
66 // for any CachedElement not in the manual clear-list below.
67 // This bit users of the SBTL surface sampler in particular,
68 // which reuses one AbstractState across thousands of (T, p)
69 // cells and was getting a constant speed_sound surface as a
70 // result.
71 cache.clear();
72
73 this->_T = -_HUGE;
74 this->_p = -_HUGE;
75 this->_Q = -_HUGE;
76 this->_rhomass.clear();
77 this->_hmass.clear();
78 this->_smass.clear();
79 this->_reverse.clear();
80 this->_phase = iphase_not_imposed;
81 return true;
82 };
83
84 // Set phase based on cached values in _p & _T or, if reverse, from (_p,_smass) or (_p,_hmass)
85 void set_phase() {
86 double epsilon = 3.3e-5; // IAPWS-IF97 RMS saturated pressure inconsistency
87 if (!_reverse) {
88 if ((std::abs(_T - IF97::Tcrit) < epsilon / 10.0) && // RMS temperature inconsistency ~ epsilon/10
89 (std::abs(_p - IF97::Pcrit) < epsilon)) { // within epsilon of [Tcrit,Pcrit]
90 _phase = iphase_critical_point; // at critical point
91 } else if (_T > IF97::Tcrit) { // to the right of the critical point
92 if (_p > IF97::Pcrit) { // above the critical point pressure
93 _phase = iphase_supercritical;
94 } else { // below the critical point pressure
95 _phase = iphase_supercritical_gas;
96 }
97 } else { // to the left of the critical point
98 if (_p > IF97::Pcrit) { // above the critical point pressure
99 _phase = iphase_supercritical_liquid;
100 } else { // at or below critical point pressure
101 double psat = IF97::psat97(_T);
102 if (_p > psat * (1.0 + epsilon)) { // above the saturation curve
103 _phase = iphase_liquid;
104 } else if (_p < psat * (1.0 - epsilon)) { // below the saturation curve
105 _phase = iphase_gas;
106 } else // exactly on saturation curve (within 1e-4 %)
107 _phase = iphase_twophase;
108 }
109 }
110 } else { // Backwards: Determine phase from _p & _smass or _hmass
111 int IF97Region;
112 if (_smass) { // Get IF97 Region
113 IF97Region = IF97::BackwardRegion(_p, _smass, IF97_SMASS); // using p, s
114 } else {
115 IF97Region = IF97::BackwardRegion(_p, _hmass, IF97_HMASS); // using p, h
116 }
117 switch (IF97Region) { // Convert IF97 Region to CP iPhase
118 case 1: // IF97::REGION1
119 if (_p <= IF97::Pcrit) {
120 _phase = iphase_liquid;
121 } else {
122 _phase = iphase_supercritical_liquid;
123 };
124 break;
125 case 2: // IF97::REGION2
126 if (_T <= IF97::Tcrit) {
127 _phase = iphase_gas;
128 } else {
129 _phase = iphase_supercritical_gas;
130 };
131 break;
132 case 3: // IF97::REGION3
133 if (_T < IF97::Tsat97(_p)) {
134 if (_p <= IF97::Pcrit) {
135 _phase = iphase_liquid;
136 } else {
137 _phase = iphase_supercritical_liquid;
138 };
139 } else {
140 if (_T <= IF97::Tcrit) {
141 _phase = iphase_gas;
142 } else {
143 _phase = iphase_supercritical_gas;
144 }
145 };
146 break;
147 case 4: // IF97::REGION4 (Saturation
148 _phase = iphase_twophase; // already handled but here for completeness
149 break;
150 case 5: // IF97::REGION5 or O.B.
151 default:
152 throw CoolProp::OutOfRangeError("Outside of IF97 Reverse Function Bounds");
153 break;
154 };
155 }
156 };
157
159
167 void update(CoolProp::input_pairs input_pair, double value1, double value2) override {
168
169 double H, S, hLmass, hVmass, sLmass, sVmass;
170
171 clear(); //clear the few cached values we are using
172
173 switch (input_pair) {
174 case PT_INPUTS:
175 _p = value1;
176 _T = value2;
177 _Q = -1;
178 set_phase();
179 // Two-Phase Check, with PT Inputs. set_phase()'s
180 // ε=3.3e-5 band around the saturation curve catches
181 // genuinely-ambiguous (p, T)-on-the-dome usage where
182 // the caller should be using PQ / QT inputs instead.
183 // But sampling-driven callers (SBTL surface build,
184 // Newton refinement of a forward h(T, p) = h_target
185 // iteration) legitimately land in this band through
186 // round-off and want a single-phase evaluation on
187 // the side they tell us. Honour an explicit
188 // specify_phase() hint here; throw otherwise.
189 if (_phase == iphase_twophase) {
190 if (imposed_phase_index == iphase_liquid || imposed_phase_index == iphase_supercritical_liquid
191 || imposed_phase_index == iphase_gas || imposed_phase_index == iphase_supercritical_gas) {
192 _phase = imposed_phase_index;
193 } else {
194 throw ValueError(format("Saturation pressure [%g Pa] corresponding to T [%g K] is within 3.3e-3 %% of given p [%Lg Pa]",
195 IF97::psat97(_T), _T, _p));
196 }
197 }
198 break;
199 case PQ_INPUTS:
200 _p = value1;
201 _Q = value2;
202 if ((_Q < 0) || (_Q > 1)) throw CoolProp::OutOfRangeError("Input vapor quality [Q] must be between 0 and 1");
203 _T = IF97::Tsat97(_p); // ...will throw exception if _P not on saturation curve
204 _phase = iphase_twophase;
205 break;
206 case QT_INPUTS:
207 _Q = value1;
208 _T = value2;
209 if ((_Q < 0) || (_Q > 1)) throw CoolProp::OutOfRangeError("Input vapor quality [Q] must be between 0 and 1");
210 _p = IF97::psat97(_T); // ...will throw exception if _P not on saturation curve
211 _phase = iphase_twophase;
212 break;
213 case HmolarP_INPUTS:
214 // IF97 is mass based so convert hmolar input to hmass
215 _hmass = value1 / molar_mass(); // Convert to mass basis : (J/mol) / (kg/mol) = J/kg
216 _p = value2;
217 // Fall thru to mass basis inputs
218 case HmassP_INPUTS:
219 if (!(_hmass)) _hmass = value1; // don't set if already set above
220 _p = value2;
221 _reverse = 1.0; // only ever set for HP and SP Inputs
222 _T = IF97::T_phmass(_p, _hmass);
223 _Q = -1; // Default if not set in 2-Phase Region
224 // ...if in the vapor dome (Region 4), calculate Quality...
225 if (IF97::BackwardRegion(_p, _hmass, IF97_HMASS) == 4) {
226 H = _hmass;
227 hVmass = IF97::hvap_p(_p);
228 hLmass = IF97::hliq_p(_p);
229 _Q = std::min(1.0, std::max(0.0, (H - hLmass) / (hVmass - hLmass))); //bound between 0 and 1
230 _phase = iphase_twophase;
231 } else {
232 set_phase();
233 };
234 break;
235 case PSmolar_INPUTS:
236 // IF97 is mass based so convert smolar input to smass
237 _p = value1;
238 _smass = value2 / molar_mass(); // Convert to mass basis : (J/mol-K) / (kg/mol) = J/kg-K
239 // Fall thru to mass basis inputs
240 case PSmass_INPUTS:
241 _p = value1;
242 if (!(_smass)) _smass = value2;
243 _reverse = 1.0;
244 _T = IF97::T_psmass(_p, _smass);
245 _Q = -1;
246 // ...if in the vapor dome (Region 4), calculate Quality...
247 if (IF97::BackwardRegion(_p, _smass, IF97_SMASS) == 4) {
248 S = _smass;
249 sVmass = IF97::svap_p(_p);
250 sLmass = IF97::sliq_p(_p);
251 _Q = std::min(1.0, std::max(0.0, (S - sLmass) / (sVmass - sLmass))); //bound between 0 and 1
252 _phase = iphase_twophase;
253 } else {
254 set_phase();
255 };
256 break;
257 case HmolarSmolar_INPUTS:
258 // IF97 is mass based so convert smolar input to smass
259 _hmass = value1 / molar_mass(); // Convert to mass basis : (J/mol) / (kg/mol) = J/kg
260 _smass = value2 / molar_mass(); // Convert to mass basis : (J/mol-K) / (kg/mol) = J/kg-K
261 // Fall thru to mass basis inputs
262 case HmassSmass_INPUTS:
263 _hmass = value1;
264 _smass = value2;
265 _p = IF97::p_hsmass(_hmass, _smass);
266 _T = IF97::T_phmass(_p, _hmass);
267 // ...if in the vapor dome (Region 4), calculate Quality...
268 if (IF97::BackwardRegion(_p, _hmass, IF97_HMASS) == 4) {
269 H = _hmass;
270 hVmass = IF97::hvap_p(_p);
271 hLmass = IF97::hliq_p(_p);
272 _Q = std::min(1.0, std::max(0.0, (H - hLmass) / (hVmass - hLmass))); //bount between 0 and 1
273 _phase = iphase_twophase;
274 } else {
275 _Q = -1;
276 set_phase();
277 };
278 break;
279 default:
280 throw ValueError("This pair of inputs is not yet supported");
281 }
282 };
283
284 /* We have to override some of the functions from the AbstractState.
285 * IF97 is only mass-based and does not support conversion
286 * from mass- to molar-specific quantities.
287 */
288 // ************************************************************************* //
289 // Basic Thermodynamic Functions //
290 // ************************************************************************* //
291 //
292 double calc_SatLiquid(parameters iCalc) {
293 switch (iCalc) {
294 case iDmass:
295 return IF97::rholiq_p(_p);
296 break;
297 case iHmass:
298 return IF97::hliq_p(_p);
299 break;
300 case iSmass:
301 return IF97::sliq_p(_p);
302 break;
303 case iCpmass:
304 return IF97::cpliq_p(_p);
305 break;
306 case iCvmass:
307 return IF97::cvliq_p(_p);
308 break;
309 case iUmass:
310 return IF97::uliq_p(_p);
311 break;
312 case ispeed_sound:
313 return IF97::speed_soundliq_p(_p);
314 break;
315 case iviscosity:
316 return IF97::viscliq_p(_p);
317 break;
318 case iconductivity:
319 return IF97::tcondliq_p(_p);
320 break;
321 case isurface_tension:
322 return IF97::sigma97(_T);
323 break;
324 case iPrandtl:
325 return IF97::prandtlliq_p(_p);
326 break;
327 default:
328 return -_HUGE;
329 };
330 }
331 double calc_SatVapor(parameters iCalc) {
332 switch (iCalc) {
333 case iDmass:
334 return IF97::rhovap_p(_p);
335 break;
336 case iHmass:
337 return IF97::hvap_p(_p);
338 break;
339 case iSmass:
340 return IF97::svap_p(_p);
341 break;
342 case iCpmass:
343 return IF97::cpvap_p(_p);
344 break;
345 case iCvmass:
346 return IF97::cvvap_p(_p);
347 break;
348 case iUmass:
349 return IF97::uvap_p(_p);
350 break;
351 case ispeed_sound:
352 return IF97::speed_soundvap_p(_p);
353 break;
354 case iviscosity:
355 return IF97::viscvap_p(_p);
356 break;
357 case iconductivity:
358 return IF97::tcondvap_p(_p);
359 break;
360 case isurface_tension:
361 return IF97::sigma97(_T);
362 break;
363 case iPrandtl:
364 return IF97::prandtlvap_p(_p);
365 break;
366 default:
367 return -_HUGE;
368 };
369 }
370 double calc_Flash(parameters iCalc) {
371 switch (_phase) {
372 case iphase_twophase: // In saturation envelope
373 if (std::abs(_Q) < 1e-10) {
374 return calc_SatLiquid(iCalc); // bubble point (Q == 0) on Sat. Liquid curve
375 } else if (std::abs(_Q - 1) < 1e-10) {
376 return calc_SatVapor(iCalc); // dew point (Q == 1) on Sat. Vapor curve
377 } else { // else "inside" bubble ( 0 < Q < 1 )
378 switch (iCalc) {
379 case iDmass:
380 // Density is an inverse phase weighted property, since it's the inverse of specific volume
381 return 1.0 / (_Q / calc_SatVapor(iDmass) + (1.0 - _Q) / calc_SatLiquid(iDmass));
382 break;
383 case iCpmass:
384 throw NotImplementedError(format("Isobaric Specific Heat not valid in two phase region"));
385 break;
386 case iCvmass:
387 throw NotImplementedError(format("Isochoric Specific Heat not valid in two phase region"));
388 break;
389 case ispeed_sound:
390 throw NotImplementedError(format("Speed of Sound not valid in two phase region"));
391 break;
392 case iviscosity:
393 throw NotImplementedError(format("Viscosity not valid in two phase region"));
394 break;
395 case iconductivity:
396 throw NotImplementedError(format("Conductivity not valid in two phase region"));
397 break;
398 case isurface_tension:
399 return IF97::sigma97(_T);
400 break; // Surface Tension is not a phase-weighted property
401 case iPrandtl:
402 throw NotImplementedError(format("Prandtl number is not valid in two phase region"));
403 break;
404 default:
405 return _Q * calc_SatVapor(iCalc) + (1.0 - _Q) * calc_SatLiquid(iCalc); // Phase weighted combination
406 };
407 }
408 break;
409 default: // Outside saturation envelope (iphase_not_imposed), let IF97 determine phase/region
410 switch (iCalc) {
411 case iDmass:
412 return IF97::rhomass_Tp(_T, _p);
413 break;
414 case iHmass:
415 return IF97::hmass_Tp(_T, _p);
416 break;
417 case iSmass:
418 return IF97::smass_Tp(_T, _p);
419 break;
420 case iCpmass:
421 return IF97::cpmass_Tp(_T, _p);
422 break;
423 case iCvmass:
424 return IF97::cvmass_Tp(_T, _p);
425 break;
426 case iUmass:
427 return IF97::umass_Tp(_T, _p);
428 break;
429 case ispeed_sound:
430 return IF97::speed_sound_Tp(_T, _p);
431 break;
432 case iviscosity:
433 return IF97::visc_Tp(_T, _p);
434 break;
435 case iconductivity:
436 return IF97::tcond_Tp(_T, _p);
437 break;
438 case isurface_tension:
439 throw NotImplementedError(format("Surface Tension is only valid within the two phase region; Try PQ or QT inputs"));
440 break;
441 case iPrandtl:
442 return IF97::prandtl_Tp(_T, _p);
443 break;
444 default:
445 throw NotImplementedError(format("Output variable not implemented in IF97 Backend"));
446 };
447 }
448 }
450 double rhomass() {
451 return calc_rhomass();
452 };
453 double calc_rhomass() override {
454 return calc_Flash(iDmass);
455 };
457 double rhomolar() {
458 return calc_rhomolar();
459 };
460 double calc_rhomolar() override {
461 return rhomass() / molar_mass();
462 };
463
465 double hmass() {
466 return calc_hmass();
467 };
468 double calc_hmass() override {
469 if (_reverse && _smass)
470 return IF97::hmass_psmass(_p, _smass); // Special IF97 function for handling reverse h(p,s) evaluation
471 else
472 return calc_Flash(iHmass);
473 };
475 double hmolar() {
476 return calc_hmolar();
477 };
478 double calc_hmolar() override {
479 return hmass() * molar_mass();
480 };
481
483 double smass() {
484 return calc_smass();
485 };
486 double calc_smass() override {
487 if (_reverse && _hmass)
488 return IF97::smass_phmass(_p, _hmass); // Special IF97 function for handling reverse s(p,h) evaluation
489 else
490 return calc_Flash(iSmass);
491 };
493 double smolar() {
494 return calc_smolar();
495 };
496 double calc_smolar() override {
497 return smass() * molar_mass();
498 };
499
501 double umass() {
502 return calc_umass();
503 };
504 double calc_umass() override {
505 return calc_Flash(iUmass);
506 };
508 double umolar() {
509 return calc_umolar();
510 };
511 double calc_umolar() override {
512 return umass() * molar_mass();
513 };
514
516 double cpmass() {
517 return calc_cpmass();
518 };
519 double calc_cpmass() override {
520 return calc_Flash(iCpmass);
521 };
523 double cpmolar() {
524 return calc_cpmolar();
525 };
526 double calc_cpmolar() override {
527 return cpmass() * molar_mass();
528 };
529
531 double cvmass() {
532 return calc_cvmass();
533 };
534 double calc_cvmass() override {
535 return calc_Flash(iCvmass);
536 };
538 double cvmolar() {
539 return calc_cvmolar();
540 };
541 double calc_cvmolar() override {
542 return cvmass() * molar_mass();
543 };
544
546 double speed_sound() {
547 return calc_speed_sound();
548 };
549 double calc_speed_sound() override {
550 return calc_Flash(ispeed_sound);
551 };
552
553 // Return the phase
554 phases calc_phase() override {
555 return _phase;
556 };
557
558 // Phase-hint plumbing. AbstractState's default calc_specify_phase
559 // throws NotImplementedError; override here so callers can use
560 // specify_phase() to disambiguate (p, T) inputs that land within
561 // IF97's 3.3e-5 ε-band of the saturation curve. The hint is read
562 // by update(PT_INPUTS) above.
563 void calc_specify_phase(phases phase) override {
564 imposed_phase_index = phase;
565 }
566 void calc_unspecify_phase() override {
567 imposed_phase_index = iphase_not_imposed;
568 }
569
570 //
571 // ************************************************************************* //
572 // Trivial Functions //
573 // ************************************************************************* //
574 //
576 double calc_Ttriple() override {
577 return IF97::get_Ttrip();
578 };
580 double calc_p_triple() override {
581 return IF97::get_ptrip();
582 };
584 double calc_T_critical() override {
585 return IF97::get_Tcrit();
586 };
588 double calc_p_critical() override {
589 return IF97::get_pcrit();
590 };
593 double calc_gas_constant() override {
594 return IF97::get_Rgas() * molar_mass();
595 };
597 double calc_molar_mass() override {
598 return IF97::get_MW();
599 };
601 double calc_acentric_factor() override {
602 return IF97::get_Acentric();
603 };
605 // TODO: May want to adjust this based on _T, since Region 5
606 // is limited to 50 MPa, instead of 100 MPa elsewhere.
607 double calc_pmax() override {
608 return IF97::get_Pmax();
609 };
612 double calc_Tmax() override {
613 return IF97::get_Tmax();
614 };
616 double calc_Tmin() override {
617 return IF97::get_Tmin();
618 };
621 double rhomolar_critical() {
622 return calc_rhomass_critical() / molar_mass();
623 }
624 double rhomass_critical() {
625 return calc_rhomass_critical();
626 }
627 // Overwrite the virtual calc_ functions for density
628 double calc_rhomolar_critical() override {
629 return rhomass_critical() / molar_mass();
630 };
631 double calc_rhomass_critical() override {
632 return IF97::get_rhocrit();
633 };
634 //
635 // ************************************************************************* //
636 // Saturation Functions //
637 // ************************************************************************* //
638 //
639 double calc_pressure() override {
640 return _p;
641 };
642 //
643 // ************************************************************************* //
644 // Transport Property Functions //
645 // ************************************************************************* //
646 //
647 // Return viscosity in [Pa-s]
648 double viscosity() {
649 return calc_viscosity();
650 };
651 double calc_viscosity() override {
652 return calc_Flash(iviscosity);
653 };
654 // Return thermal conductivity in [W/m-K]
655 double conductivity() {
656 return calc_conductivity();
657 };
658 double calc_conductivity() override {
659 return calc_Flash(iconductivity);
660 };
661 // Return surface tension in [N/m]
662 double surface_tension() {
663 return calc_surface_tension();
664 };
665 double calc_surface_tension() override {
666 return calc_Flash(isurface_tension);
667 };
668 // Return Prandtl number (mu*Cp/k) [dimensionless]
669 double Prandtl() {
670 return calc_Flash(iPrandtl);
671 };
672
677 void fast_evaluate(CoolProp::input_pairs input_pair, const double* val1, const double* val2, std::size_t N_inputs,
678 const CoolProp::parameters* outputs, std::size_t N_outputs, double* out_buffer, std::size_t out_buffer_size, int* status_flags,
679 std::size_t status_flags_size, CoolProp::phases imposed_phase = CoolProp::iphase_not_imposed) override {
680 if (N_outputs != 0 && N_inputs > std::numeric_limits<std::size_t>::max() / N_outputs) {
681 throw ValueError(format("fast_evaluate: N_inputs * N_outputs would overflow size_t (N_inputs=%zu, N_outputs=%zu)", N_inputs, N_outputs));
682 }
683 const std::size_t required_out = N_inputs * N_outputs;
684 if (out_buffer_size < required_out) {
685 throw ValueError(format("fast_evaluate: out_buffer_size=%zu < required %zu (N_inputs * N_outputs)", out_buffer_size, required_out));
686 }
687 if (status_flags_size < N_inputs) {
688 throw ValueError(format("fast_evaluate: status_flags_size=%zu < required %zu (N_inputs)", status_flags_size, N_inputs));
689 }
690 if (N_inputs == 0) return;
691 // Null-pointer check BEFORE the N_outputs==0 early-write path: that path
692 // dereferences status_flags, so the check has to dominate it.
693 if (val1 == nullptr || val2 == nullptr || outputs == nullptr || out_buffer == nullptr || status_flags == nullptr) {
694 throw ValueError("fast_evaluate: null pointer argument");
695 }
696 if (N_outputs == 0) {
697 for (std::size_t k = 0; k < N_inputs; ++k)
698 status_flags[k] = CoolProp::fast_evaluate_ok;
699 return;
700 }
701
702 const bool is_PT = (input_pair == PT_INPUTS);
703 const bool is_HmassP = (input_pair == HmassP_INPUTS);
704 const bool is_HmolarP = (input_pair == HmolarP_INPUTS);
705 if (!is_PT && !is_HmassP && !is_HmolarP) {
706 throw ValueError(format("fast_evaluate (IF97): input_pair %s not supported (use PT_INPUTS, HmassP_INPUTS, HmolarP_INPUTS)",
707 get_input_pair_short_desc(input_pair).c_str()));
708 }
709
710 // Validate output keys up-front. IF97 covers everything its calc_Flash supports.
711 constexpr std::size_t MAX_FE_OUTPUTS = 64;
712 if (N_outputs > MAX_FE_OUTPUTS) {
713 throw ValueError(format("fast_evaluate: N_outputs=%zu exceeds compile-time limit %zu", N_outputs, MAX_FE_OUTPUTS));
714 }
715 for (std::size_t o = 0; o < N_outputs; ++o) {
716 switch (outputs[o]) {
717 case iT:
718 case iP:
719 case iDmass:
720 case iDmolar:
721 case iHmass:
722 case iHmolar:
723 case iSmass:
724 case iSmolar:
725 case iUmass:
726 case iUmolar:
727 case iCpmass:
728 case iCpmolar:
729 case iCvmass:
730 case iCvmolar:
731 case ispeed_sound:
732 case iviscosity:
733 case iconductivity:
734 case iPrandtl:
735 case iQ:
736 break;
737 default:
738 throw ValueError(
739 format("fast_evaluate (IF97): output[%zu]=%s not supported", o, get_parameter_information(outputs[o], "short").c_str()));
740 }
741 }
742
743 const double M = calc_molar_mass(); // kg/mol, constant for water
744 const double NaN = std::numeric_limits<double>::quiet_NaN();
745
746 // imposed_phase: IF97 derives phase from (T,p) — we honour the hint by
747 // skipping the two-phase check when caller asserts single-phase.
748 const bool skip_twophase_check = (imposed_phase != iphase_not_imposed && imposed_phase != iphase_twophase && imposed_phase != iphase_unknown);
749
750 auto fill_nan_row = [&](std::size_t k) {
751 for (std::size_t o = 0; o < N_outputs; ++o) {
752 out_buffer[k * N_outputs + o] = NaN;
753 }
754 };
755
756 for (std::size_t k = 0; k < N_inputs; ++k) {
757 const double v1 = val1[k];
758 const double v2 = val2[k];
759
760 double T_k, p_k;
761 const double hmass_k = is_HmassP ? v1 : (is_HmolarP ? (v1 / M) : 0.0);
762 bool dome_hit = false;
763 double T_sat_k = 0.0, Q_k = 0.0, hL_k = 0.0, hV_k = 0.0;
764 if (is_PT) {
765 p_k = v1;
766 T_k = v2;
767 } else {
768 // (h, p) input — convert h to mass basis and use IF97's
769 // backward T_phmass to recover T. Both branches must produce
770 // a clean (T, p) for the forward evaluators below.
771 p_k = v2;
772 bool inverse_ok = false;
773 try {
774 T_k = IF97::T_phmass(p_k, hmass_k);
775 inverse_ok = true;
776 } catch (const std::exception&) { // NOLINT(bugprone-empty-catch)
777 // T_phmass throws for out-of-range inputs — fall
778 // through to the dome probe + OOB classification.
779 }
780 // T_phmass returns Tsat for in-dome (h, p) inputs without
781 // throwing, so the dome check has to run independently of
782 // whether the inverse "succeeded". Probing hL / hV at p
783 // and bracketing on h is the cleanest detector.
784 if (p_k > 0 && p_k <= IF97::Pcrit) {
785 try {
786 const double hL_probe = IF97::hliq_p(p_k);
787 const double hV_probe = IF97::hvap_p(p_k);
788 if (hmass_k >= hL_probe && hmass_k <= hV_probe) {
789 T_sat_k = IF97::Tsat97(p_k);
790 hL_k = hL_probe;
791 hV_k = hV_probe;
792 Q_k = (hmass_k - hL_k) / (hV_k - hL_k);
793 dome_hit = true;
794 T_k = T_sat_k;
795 inverse_ok = true;
796 }
797 } catch (const std::exception&) { // NOLINT(bugprone-empty-catch)
798 // Above pcrit or other IF97 envelope failure;
799 // fall through to OOB.
800 }
801 }
802 if (!inverse_ok) {
803 status_flags[k] = CoolProp::fast_evaluate_out_of_range;
804 fill_nan_row(k);
805 continue;
806 }
807 }
808
809 // Range check before evaluation. IF97 has explicit T/p limits; we
810 // reject points outside and let in-range points go to the kernel,
811 // which itself will throw on Region-5 etc. (caught below).
812 if (!(p_k > 0 && T_k > 0 && p_k <= IF97::get_Pmax() && T_k >= IF97::get_Tmin() && T_k <= IF97::get_Tmax())) {
813 status_flags[k] = CoolProp::fast_evaluate_out_of_range;
814 fill_nan_row(k);
815 continue;
816 }
817
818 // Two-phase rejection: PT_INPUTS landing exactly on the saturation
819 // curve is ambiguous; caller can use imposed_phase or QT/PQ inputs
820 // (which we don't support here yet). For (h,p) inputs the backward
821 // solver places us in a specific region so this is less common.
822 if (!skip_twophase_check && is_PT && T_k <= IF97::Tcrit) {
823 const double psat = IF97::psat97(T_k);
824 const double eps = 3.3e-5;
825 if (std::abs(p_k - psat) <= psat * eps) {
826 status_flags[k] = CoolProp::fast_evaluate_two_phase_disallowed;
827 fill_nan_row(k);
828 continue;
829 }
830 }
831
832 bool eval_failed = false;
833 // Two-phase Q-blend cache (lazy fill — only touched if a
834 // property that needs sat-line endpoints is requested while
835 // dome_hit is true).
836 double rhoL_k = 0.0, rhoV_k = 0.0;
837 double sL_k = 0.0, sV_k = 0.0;
838 double uL_k = 0.0, uV_k = 0.0;
839 bool have_rho = false, have_s = false, have_u = false;
840 auto ensure_rho_sat = [&] {
841 if (have_rho) return;
842 rhoL_k = IF97::rholiq_p(p_k);
843 rhoV_k = IF97::rhovap_p(p_k);
844 have_rho = true;
845 };
846 auto ensure_s_sat = [&] {
847 if (have_s) return;
848 sL_k = IF97::sliq_p(p_k);
849 sV_k = IF97::svap_p(p_k);
850 have_s = true;
851 };
852 auto ensure_u_sat = [&] {
853 if (have_u) return;
854 uL_k = IF97::uliq_p(p_k);
855 uV_k = IF97::uvap_p(p_k);
856 have_u = true;
857 };
858 for (std::size_t o = 0; o < N_outputs; ++o) {
859 const parameters out_key = outputs[o];
860 try {
861 double val;
862 if (dome_hit) {
863 // Q-weighted blend for two-phase points. Specific
864 // volume blends linearly (density is its reciprocal);
865 // h/s/u blend linearly directly. Speed of sound,
866 // cp, cv, viscosity, conductivity have no
867 // equilibrium bulk value in the dome — set NaN.
868 switch (out_key) {
869 case iT:
870 val = T_sat_k;
871 break;
872 case iP:
873 val = p_k;
874 break;
875 case iHmass:
876 val = hmass_k;
877 break;
878 case iHmolar:
879 val = hmass_k * M;
880 break;
881 case iDmass:
882 case iDmolar: {
883 ensure_rho_sat();
884 const double vL = 1.0 / rhoL_k;
885 const double vV = 1.0 / rhoV_k;
886 const double rho_mass = 1.0 / ((1.0 - Q_k) * vL + Q_k * vV);
887 val = (out_key == iDmass) ? rho_mass : rho_mass / M;
888 break;
889 }
890 case iSmass:
891 case iSmolar: {
892 ensure_s_sat();
893 const double s_mass = (1.0 - Q_k) * sL_k + Q_k * sV_k;
894 val = (out_key == iSmass) ? s_mass : s_mass * M;
895 break;
896 }
897 case iUmass:
898 case iUmolar: {
899 ensure_u_sat();
900 const double u_mass = (1.0 - Q_k) * uL_k + Q_k * uV_k;
901 val = (out_key == iUmass) ? u_mass : u_mass * M;
902 break;
903 }
904 case iQ:
905 val = Q_k;
906 break;
907 case iCpmass:
908 case iCpmolar:
909 case iCvmass:
910 case iCvmolar:
911 case ispeed_sound:
912 case iviscosity:
913 case iconductivity:
914 case iPrandtl:
915 val = NaN;
916 break;
917 default:
918 val = NaN;
919 eval_failed = true;
920 break;
921 }
922 } else
923 switch (out_key) {
924 case iT:
925 val = T_k;
926 break;
927 case iP:
928 val = p_k;
929 break;
930 case iDmass:
931 val = IF97::rhomass_Tp(T_k, p_k);
932 break;
933 case iDmolar:
934 val = IF97::rhomass_Tp(T_k, p_k) / M;
935 break;
936 case iHmass:
937 val = IF97::hmass_Tp(T_k, p_k);
938 break;
939 case iHmolar:
940 val = IF97::hmass_Tp(T_k, p_k) * M;
941 break;
942 case iSmass:
943 val = IF97::smass_Tp(T_k, p_k);
944 break;
945 case iSmolar:
946 val = IF97::smass_Tp(T_k, p_k) * M;
947 break;
948 case iUmass:
949 val = IF97::umass_Tp(T_k, p_k);
950 break;
951 case iUmolar:
952 val = IF97::umass_Tp(T_k, p_k) * M;
953 break;
954 case iCpmass:
955 val = IF97::cpmass_Tp(T_k, p_k);
956 break;
957 case iCpmolar:
958 val = IF97::cpmass_Tp(T_k, p_k) * M;
959 break;
960 case iCvmass:
961 val = IF97::cvmass_Tp(T_k, p_k);
962 break;
963 case iCvmolar:
964 val = IF97::cvmass_Tp(T_k, p_k) * M;
965 break;
966 case ispeed_sound:
967 val = IF97::speed_sound_Tp(T_k, p_k);
968 break;
969 case iviscosity:
970 val = IF97::visc_Tp(T_k, p_k);
971 break;
972 case iconductivity:
973 val = IF97::tcond_Tp(T_k, p_k);
974 break;
975 case iPrandtl:
976 val = IF97::prandtl_Tp(T_k, p_k);
977 break;
978 case iQ:
979 // Single-phase by definition here — the
980 // dome-detection branch above would have
981 // routed two-phase points to the
982 // Q-blend path. Mirror update()'s
983 // -1 sentinel for single-phase Q.
984 val = -1.0;
985 break;
986 default:
987 val = NaN;
988 eval_failed = true;
989 break;
990 }
991 out_buffer[k * N_outputs + o] = val;
992 } catch (const std::exception&) {
993 eval_failed = true;
994 out_buffer[k * N_outputs + o] = NaN;
995 }
996 }
997 // API contract: when status_flags[k] is non-zero, the entire row is
998 // NaN. The per-output catch above only NaNs the specific failing
999 // output; finish the job here so callers don't see partial rows.
1000 if (eval_failed) {
1001 fill_nan_row(k);
1002 status_flags[k] = CoolProp::fast_evaluate_internal_error;
1003 } else {
1004 status_flags[k] = CoolProp::fast_evaluate_ok;
1005 }
1006 }
1007 // IF97 is stateless externally; clear() resets internal cached values to
1008 // be consistent with the TabularBackend contract (state untouched).
1009 clear();
1010 };
1011};
1012
1013} /* namespace CoolProp */
1014#endif /* IF97BACKEND_H_ */
Generated by doxygen 1.9.4